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ABSTRACT 


As part of the study being made at the N.A.C.A. Cleveland 
laboratory on the flutter of compressor and turbine blades, a 
theoretical analysis was made of the effect of aerodynamic hyster- 
esis on stalling flutter. The assumption was made that the 
absolute magnitude of the oscillatory aerodynamic forces and 
moments are the same at stall as at zero angle of attack but that 
the vector magnitudes of these forces and moments are changed, 
this change being caused by the lag of aerodynamic damping and 
restoring forces behind the velocities and displacements at stall, 
The decrease of critical 
The results 


thus giving rise to a hysteresis effect. 
flutter speed at stall was thus theoretically shown. 
were applied to a given airfoil, and correlation of the experimental 
and theoretical results was found possible by assuming that the 
angle of aerodynamic lag varies as the slope of the static-lift 
curve. The aerodynamic lag was shown to cause the effective 
torsional damping to decrease, thereby explaining the low values 
of torsional aerodynamic damping obtained at stall. 


INTRODUCTION 


~~ THE LARGE INCREASE in the use of axial-flow 
compressors and turbines, an investigation of the 
vibrations of the blades of these units has become of in- 
creasing importance. A significant part of such a study 
involves the type of vibration that is self-maintained by 
the continual absorption of energy from the air stream 
namely, the flutter of turbine and compressor blades. 
The term ‘‘flutter,’’ as it is ordinarily used, refers to 
classical flutter, a self-sustained oscillation due to the 
coupling of inertia forces, elastic forces, damping forces, 
and dynamic aerodynamic forces. This type of flutter 
usually occurs on airplane wings at low angles of attack 
When the velocity reaches a certain value determined 
by the wing design and is called the critical flutter 
speed. The theory for this type of flutter has been 
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developed by many investigators, and the oscillatory 
aerodynamic forces and moments have been derived.'~* 
In all these derivations, the airfoil is assumed to be at 
zero angle of attack. Close agreement has been found 
between theory and experiment. 


A radically different type of flutter called ‘‘stalling 
flutter,’’ however, may occur on airfoils at high angles 
of attack, such as blades near the stall point. This 
flutter occurs at much lower critical flutter speeds 
than would be expected from classical-flutter theory. 
Two reasons exist for believing that this type of flutter 
would be more likely to occur in compressor and tur 
bine blades than classical flutter: First, compressor and 
turbine blades are so stiff compared with ordinary air 
plane wings that critical speeds obtained by classical 
calculations are usually above the reasonable operating 
range, at least for subsonic airfoils. This conclusion, 
however, cannot be stated with certainty because the 
effect of cascading on the oscillatory aerodynamic forces 
is as yet unknown. Second, the operation of com- 
pressors and turbines at high loads near the stall point 
makes the occurrence of stalling flutter more likely. 
Stalling flutter has been shown to occur on compressor 
cascades in bench tests.‘ 


Experimental investigations carried out at the N.A. 
C.A. Cleveland laboratory on compressors under oper 
ating conditions have shown that in some cases the vi- 
bratory stresses increase greatly with angle of attack. 
Inasmuch as these vibrations were apparently aerody- 
namically excited, it was decided that a basic under- 
standing of the phenomenon of flutter at high angles 
of attack would be needed before the problem could 
be successfully solved. The investigation of flutter 
at stall, part of which is described in the present paper, 
was therefore undertaken. 
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Angle of attack, a, deg 
Fic. 1. Typical lift curve showing hysteresis loop. 

As yet, little is known about the phenomenon of 
stalling flutter. It has been observed most frequently 
in connection with stalled propellers.6 Experimental 
work has shown three possible causes.°-* The first 
cause is termed ‘‘static instability’’® and is associated 
with the fact that the slope of the lift curve decreases 
and becomes negative near the stall point. This change 
of slope can cause an instability, which has been ex- 
plained by Den Hartog® as follows: As the airfoil 
moves upward, the resultant motion decreases the ef- 
fective angle of attack. If the slope of the lift curve 
is negative, the lift force will increase. As the airfoil 
moves downward, the effective angle of attack in- 
creases and the lift force decreases. The vibrations 
therefore continue to build up. This explanation in 
terms of the static lift curve is, of course, not strictly 
correct for conditions of flutter because the dynamic- 
lift curve will be different from the static-lift curve. 


The second possible cause for stalling flutter shown 
to exist® was excitation of the airfoil by a system of 
K4rm4n vortices shed in the wake. Karman vortices 
produce an excitation with frequency directly propor- 
tional to the air velocity. At certain discrete velocities, 
the K4rm4n vortex frequency will coincide with one of 
the natural frequencies of the air foil. When this coin- 
cidence of frequencies occurs, resonant vibrations of the 
airfoil can be built up. 


The third and apparently most common cause was 
first investigated by Stiider,’ who carried out the first 
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important experimental work on the subject. Stiider 
found that the critical speed of his models dropped 
considerably near the stall point and that the phase 
difference between the bending and torsional oscillations 
decreased to zero so that the torsional and bending mo. 
tions were in phase. He also found that torsional one. 
degree-of-freedom flutter was possible, whereas it js 
known that one-degree-of-freedom classical flutter can- 
Stiider explained the phenomenon occurring 
As the airfoil approaches the stall 
As the airfoil os. 


not occur. 
at stall as follows: 
angle, separation of flow occurs. 
cillates about the stall point, the separation of flow js 
delayed until the position of maximum amplitude js 
reached. This position will be at an angle of attack 
greater than that for the stationary airfoil. On the 
return movement, re-establishment of smooth flow is 
delayed to an angle of attack below that at which the 
stationary airfoil would stall. At the stall angle, there- 
fore, the air forces acting depend on the direction of 
motion, giving rise to a hysteresis loop at stall, as shown 
in Fig. 1. Stiider came to the conclusion that this 
hysteresis effect, enabling energy to be absorbed 
from the air stream, is the cause of stalling flut- 
ter. 

Stiider’s work was verified and amplified by Victory,’ 
who used it as a basis for finding a link between stalling 
flutter and classical theory. The torsional aerodynamic 
damping was experimentally shown to decrease as the 
angle of attack increases, and it may become negative 
at stall. It was then shown that the critical flutter 
speed at stall could be calculated with reasonable ac- 
curacy by classical theory if the usually accepted classi- 
cal value of the torsional aerodynamic damping is re- 
placed by an experimental value that varies with angle 
of attack, frequency parameter, Reynolds Number, 
and position of elastic axis. 

A similar hysteresis phenomenon was investigated 
by Smilg in connection with aileron flutter at transonic 
speeds. There the flow separation is not due to stall 
but rather to compressibility effects in the transonic 
range. The effects, however, are similar. The pro- 
cedure used was to assume that because of the flow 
separation the damping and restoring aerody- 
namic forces and moments lag the velocities and 
displacements, respectively. This procedure  sug- 
gested a similar approach to the problem of stalling flut- 
ter. 

The mathematical theory of the aerodynamic hystere- 
sis effect making use of the concept of an aerodynamic 
angle of lag developed at the N.A.C.A. Cleveland lab- 
oratory and the results obtained on a given airfoil are 
presented herein. 


SYMBOLS 
a = coordinate of elastic axis measured from mid- 
chord, positive towards trailing edge ™ 

units of half chord 





~_ cae a 


-_ 


time 


Crit 


give 
If t] 
and 
tion 
tion 


Stiider 
ropped 
- phase 
lations 
ng mo- 
al one. 
S it is 
er can- 
curring 
1e stall 
oil os- 
flow is 
‘ude is 
attack 
Jn the 
flow is 
ich the 
there- 
tion of 
shown 
it this 
sorbed 
x flut- 


ctory,$ 
talling 
mamic 
as the 
ative 
flutter 
ple ac- 
classi- 
iS re- 
| angle 
imber, 


igated 
nsonic 
o stall 
nsonic 
e pro- 
e flow 
erody- 
s and 

sug: 
g flut- 


ystere- 
namic 
d lab- 
oil are 


m mid- 
ge if 





FLUTTER 


a;, 01, (1, @2, = functions of reduced frequency & and angle of 
be, C2 aerodynamic lag ¢ 

ay, 4,, 20,49 = aerodynamic coefficients 

b = half chord used as reference unit length 

b,, b,, be, 68 = aerodynamic coefficients 

c = ury? 

c = stiffness of airfoil in bending per unit span 
length 

G. = lift coefficient 

"sa = stiffness of airfoil in torsion about @ per unit 
span length 

e =_— 

F = function of reduced frequency & given in ref- 
erence | 

F,, Fe = functions of reduced frequency & and angle of 
aerodynamic lag ¢ 

G = function of reduced frequency & given in ref- 
erence | 

I = moment of inertia about elastic axis per unit 
span length 

=V-1 

K = proportionality constant 

k = reduced frequency wh/v 

L = oscillatory aerodynamic lift force per unit 
span 

M = oscillatory aerodynamic moment abou: a per 
unit span 

m = mass of airfoil per unit span 

mM, = mass of surrounding air cylinder per unit span, 
wpb? 

r = location of center of gravity of airfoil meas- 
ured from a in units of half chord 

> = radius of gyration referred to a in units of half 
chord 

S = static mass unbalance of airfoil, mrb 

t = time 

v = critical flutter speed 

y = vertical displacement 

Yo = maximum bending amplitude 

a = static angle of attack 

B = phase angle between bending and torsion 

6 = angle of torsional displacement 

9 = maximum torsional amplitude 

m = ratio of mass of airfoil per unit span to mass 
of surrounding air cylinder per unit span, 
m/m, =m/mrpb? 

p = density of surrounding air 

¢ = angle of lag between displacements or veloci- 
ties and aerodynamic restoring forces or 
damping forces 

w = frequency of flutter vibration 

w, = fundamental bending frequency of airfoil 

w, = fundamental torsional frequency of airfoil 


Dots above the symbols represent derivatives with respect to 
time. 


ANALYSIS 


Critical Flutter Speed and Frequency 


The oscillatory aerodynamic lift and moment are 
given by equations XVIII and XX of reference 1. 
If the exponential relations y = yoe“ and 0 = Oe" ~ ” 
and the relation k = wb/v are substituted in these equa- 
tions in order to eliminate imaginary terms, the equa- 
tions can be written as follows: 
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= l r 
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hh’ 2 ee 


b, = a+ (2/k)[a + (1/2)]G 


R 
°c 
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b, = —(2/k) [a + (1/2)]F 

i E r ( a G E 

) = _ - _— - — - 

> gtroet nets Av a’) G| 


2 \. 2/1. 1/1 
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The second derivative, the first derivative, and the 
displacement terms in Eqs. (1) represent inertia terms, 
damping terms, and stiffness terms, respectively. The 
bending and the torsional displacement, as has been al- 
ready indicated, are given by 
y= yor" 
and 

0 = Oe’ — B) 


Eqs. (1) then become 


Zz. = —maa*| aire oa idyyoe™ - (« —_ *) p 4 
bOye* et - ® fe iGigbOoe' — B) he ay bOve" - »| 
k ' 
- ) 
M= —mau'0| bye + ibyyoe™ + (6. — *) 4 


bOpe — B) +4. ib bOoe' — B) af * bO.e** - | 

These equations represent the vector sum of the aero- 
dynamic inertia, damping, and stiffness forces and 
moments. 

The separation and re-establishment of flow as the 
airfoil oscillates about the stall point have been shown’ 
to lag behind the actual airfoil displacement. The 
oscillatory aerodynamic forces and moments can be 
resolved as shown into inertia components in phase with 
the accelerations, damping components in phase 
with the velocities, and stiffness components in phase 
with the displacements. If the flow separation and re- 
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establishment lag the motions, these aerodynamic force 
and moment components lag the velocities and the 
displacements. The essential assumption is now made 
that at stall the absolute values of the oscillating aerody- 


namic forces and moments are the same as at zero angle of 


attack, but, because of the hysteresis effect, the vector mag- 
nitudes of the aerodynamic restoring forces and moments 
and the aerodynamic damping forces and moments have 
so changed that the restoring forces and moments lag 
the displacements and the damping forces and moments 
lag the velocities by an angle ¢. 

Eqs. (2) are then modified as follows: 


L=- Maw” | ean aa 1ayyoe’ - #7) aa 


a + = ae. p 
(a = *) bbe" . + 1a,b00e i 4 y i 


a, ea ee 
bAve’' t f . | 
k 


M= - maal be + iby ye’ ~ * + 


o . rie 
(2 . ; bore! ~ ©) + sBebbye et — 8-4 


h 
. bA en “ee | ’ 
k J 


The equations of dynamic equilibrium of the system 


0 
\ 
0 


—muw*yoe' + Cyyre™ — Swhe'"- 9 — L = °t (4a) 
— Sw*yoe™ a TwO ye!“ — B) }- C Bye" — ff) __ M a 0! a 


are:} 


I 


my + Cy + S6-— L 
Sy +16 +C0-M= 


or 


By combining Eqs. (4a) and (3) and dividing the first 
equation by mw? and the second equation by m,w*d in 
order to make the coefficients dimensionless, the follow- 
ing equations are obtained: 


m & ee = j 
| («. = Si +) er + ide’ Is + 
Me Mw? 
| (« . a )ei D4 (= + iay) x | (5) 
— =— , ” 
"  k  m~b k . 


} 


where 
- ~( $19 ane? « 
dy = pc(w,/w,)? (2 cos? ¢ 1) 
dy = pc(a/w,)? sin 2¢ 
ae 
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2 t.. ro ee 
by — ee + ibe ~ © Wyo + 
m,b, 
b, I OF (5 
b, a Po +. : e' wi p) A. ; \0, 
( k mb? — m,w*b? } Cont, 
a. ee 
( + 2w,) 6" " 166 = 0 
kh 
Let mm, =u. Then 
Cy, = mw,"?; C,/mqw* = p(w?/w?) 
f= mJ: 1/a0* = p77 = « 


5S = mrb; S/mb = ur =e 
. * C, Wt w," 
C, = Iw;’; — = i = "C— 
m,w"b? @- @w- 


By substituting into Eq. (5) and factoring as shown, 


wp” = = ot 
E + u( ae \) + ide « ley + 
( ay )+ dy 4 ay) ig 
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Eqs. (6) have solutions other than @ = y = 0 if and 


only if the determinant of the coefficients of y and 6 


vanishes. 
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k k 
By multiplying through by e’* and substituting the 
Euler relation 
e* =cosg +ising 
and by separating real and imaginary parts and rear- 
ranging : 
a;(w,/w)* + di(w,/w)? + 4 = 0) 
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By solving Eq. (8) for (w,;/w)* and equating, a = —0.29 
/ b = 3.375 in. 
(“) —h, - > b,? oe $ajcy = — by = \ bo? —_ bdsCe r = 0.228 
w 2a, Pay = 161.2 
(0) (10) w, = 87.2 rad. per sec. 
w, = 80.3 rad. per sec. 
or oan ° . ° . . 
Phe radius of gyration r,, which was not given, was as 
Fy (k, g) = Fr (k, ¢) (11) sumed to equal (6.5. The results are shown in Fig. 2. 
For a given airfoil, the functions F, and F) are func- As the angle of aerodynamic lag ¢ increases, the critical 
« ~ c , . c « « . . 
if and tions only of the reduced frequency & and the angle of flutter speed t dec — If the torsional aerody mane 
and § aerodynamic lag gy. For a given angle of aerodynamic damping coefficient 4, is now calculated by means of 
; b > 2\ ea — — ;?; oe 3. aina The 
lag y, Eq. (11) can be solved for reduced frequency k Eq. ( 13), the curve shown a Fig. 3 is obtained. Phe 
by finding the intersections of the functions F, and F: torsional aerodynamic damping coefficient 4, decreases 
. 5 . . . « 2 . e . > 
plotted against k. The flutter frequency w can then with ene angle of lag ¢. From Figs. 2 and 3, 
be obtained from Eq. (10) and the critical flutter speed the critical flutter speed , can be obtained as 8 function 
| vfrom the relation of the torsional damping coefficient by. This relation is 
(7) plotted in Fig. 4 and shows the critical flutter speed v 
v = wh/k (12) decreasing as the torsional damping 5, decreases. The 
The critical flutter speeds v and the frequencies & can shape of this curve is similar to those obtained in refer 
then be plotted as functions of the angle of aerodynamic —_ for ratios of torsional to bending frequency 
lag w/w» Close to 1. 
ap ¢. 
ng the 
Aerodynamic Damping 
The torsional aerodynamic damping coeflicient 4, 
isgiven by 
1 rear- ‘ 
, _ fl (: \p a(3 . )é oe 
4 = —-¢g<— Z —_— @¢ po a 7 (13) 2 
kL2 | k\2 ! 
; + 
(8) j . ; " : , & 
and is a function of k only. After the reduced fre-  * 
© te 
quency & is obtained as a function of the angle of aero- : 
dynamic lag yg by Eq. (11), the torsional aerodynamic & 
damping 5, can be calculated as a function of the angle 3 
ot aerodynamic lag ¢. g 
RESULTS AND DISCUSSION 
- The critical flutter speed v was calculated as a func- J 
| tion of the angle of aerodynamic lag ¢ by means of Eqs. 
(11\ . . ae Angle of aerod ic lag, % 4 
(11) and (12) for wing II of reference 6. The constants ee Re : Receestaly -htdngs . 
ithe wit Fic. 2. Variation of theoretical critical flutter speed with angle 
© this wing are: of aerodynamic lag. 
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FLUTTER 


The results obtained thus far indicate that the aero- 
dynamie-lag effect can account for a drop in critical 
flutter speed v in the region of stall and also for a de- 
crease in torsional aerodynamic damping. In order 
to obtain quantitative agreement with the experimen- 
tally observed critical flutter speeds in the region of stall, 
however, it is necessary to find a relation between the 
angle of aerodynamic lag ¢ and the angle of attack a. 
Fig. 7 of reference 6 shows the variation of critical 
flutter speed v with angle of attack a for this airfoil. 
By cross-plotting this curve with Fig. 2, the angle of 
aerodynamic lag ¢ can be obtained as a function of 
angle of attack a. This relation is shown plotted as a 
solid line in Fig. 5. 

The question now presents itself of whether some 
physical basis can be found for the curve shown in Fig. 
5. A relation between angle of attack a and the angle 
of aerodynamic lag ¢ suggests itself when the lift curve 
of the airfoil is considered. As the angle of attack 
a increases, the separation of flow increases and the 
slope of the lift curve decreases. Also, as the separa- 
tion of flow increases, the angle of aerodynamic lag ¢ 
might be expected to increase. It would therefore 
seem possible that the angle of aerodynamic lag ¢ 
is related to the change of slope to the lift curve. The 
lift curve for the airfoil used is given in reference 6. 
The assumption was now made that the angle of aero 
dynamic lag yg is given by 

¢ = K[(dC,/da), ~ 4, — (dCz/da)] 
The slope of the lift curve at zero angle of attack is ap- 
proximately 27, and the constant A was calculated by 
taking one arbitrary point from the solid curve of Fig. 
5. The value of K is equal to 45/27 if the angle of aero- 
dynamic lag gisin degrees. Therefore, 


(14) 


g = (45/2r) [24 — (dC,/da)] 





This equation was plotted as the dashed curve of 
Fig. § and is seen to agree very well with the solid curve 
which is based on experimental data. 


By combining the dashed curve of Fig. 5 with the re- 
lation shown in Fig. 2, the critical flutter speed v can 
now be plotted as a function of angle of attack a. This 
function is shown as the dashed curve of Fig. 6, which 
agrees very well with the experimental curve obtained 
from Fig. 7 of reference 6. (The slope of the lift curve 
shown in reference 6, by means of which the dashed 
curve of Fig. 5 was obtained, is plotted in Fig. 7 of this 
report. ) 

The torsional aerodynamic damping }, can now also 
be plotted as a function of the angle of attack a (Fig. 
8). The torsional aerodynamic damping 5, drops 
sharply in the region of stall, the shape of the curve 
being similar to the shape of the critical flutter-speed 
curve shown in Fig. 6. 
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SUMMARY OF RESULTS 


From the calculations of critical flutter speeds made 
for a given airfoil, it was shown that the phenomenon 
of stalling flutter can at least in some cases be ex- 
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plained on the basis of an aerodynamic lag or hysteresis 
effect. The general characteristics of stalling flutter 

namely, decreases in critical flutter speed and in ef- 
fective aerodynamic torsional damping coefficient with 
increasing angle of attack—were shown to be related 
to the more fundamental concept of lagging aerody- 
namic forces and moments. Correlation between ex- 
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perimental and theoretical results was obtained by as. 
suming that the aerodynamic angle of lag varies with 
the angle of attack in a manner that can be explained 
by the variation in the slope of the static-lift curve. 
More experimental and theoretical work is necessary to 
correlate more closely the aerodynamic lag with the 
angle of attack. 
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The Optimum Proportions of a Multicell 
Box Beam in Pure Bending’ 


HERBERT BECKER? 
Edo Corporation 


SUMMARY 


Employing a method similar to that used for long, unstiffened 
circular cylinders in pure bending, a theory is developed which 
shows that there is an optimum chord for a multicell box beam of 
specific structural thickness ratio resisting a constant moment. 
The theory, utilizing the secant modulus for plastic buckling of 
the compression cover, assumes that cover buckling constitutes 
failure. This assumption is applicable to highly stressed wings 
found on transonic and supersonic aircraft. 

A simplified calculation procedure is included for rapid design 


purposes. 

The investigation reveals that optimum cross-section propor- 
tions require a web spacing equal to the box depth; any box de- 
signed to fail at a stress in the neighborhood of the yield point 
will be close to the optimum design; and the optimum stress for 
a multicell box beam in pure bending is approximately the same 
as for a long, unstiffened circular cylinder under the same load- 


ing. 
INTRODUCTION 


y HAS BEEN DEMONSTRATED that a stress-strain curve 
provides sufficient information to obtain a mini- 
mum-weight design for a long, unstiffened circular 
cylinder in pure bending.' Furthermore, the theory 
leading to the design yields a unique radius for the ap- 
plied moment and material. 

In a manner similar to that used for cylinders, it is 
possible to show that a multicell box beam under pure 
moment which fails by buckling of the compression 
cover has an optimum structural chord that is a func- 
tion of the moment, material stress-strain curve, and 
structural thickness ratio. The theory leading to this 
optimum chord assumes the secant modulus to be ap- 
plicable to plate buckling in the plastic range. There 
is little difference between the weight determined on this 
premise and the weight that would be obtained if the 
plastic modulus, described by Stowell in reference 2, for 
simply supported plates were to be used. This modulus 
differs from the secant modulus by a modifying coeffi- 
cient involving the ratio of tangent to secant modulus 
and is referred to hereafter as the “‘modified secant 
modulus.”’ 

In order to obtain a minimum weight for the box 
beam, the compression cover and its stiffening webs 
must be arranged in the lightest possible manner. (The 

Received March 3, 1949. 

* The author wishes to express his indebtedness to Leon Kaplan 
of the Edo Corporation for his helpful and stimulating sugges- 
tions and criticisms, especially with regard to the discussion of 
the theory. 

t Now, Structures Engineer, Combustion Engineering—Super- 
heater, Marine Department. 


tension cover usually will be designed to a constant 
stress and, consequently, is not a variable quantity. 
Therefore its effect will be disregarded, and all refer- 
ences to weight will imply that of compression cover 
plus webs.) It was found by Gerard that this would 
occur for a definite relationship between the number of 
bays and the value of the structural thickness ratio. 
This relationship is expressed in Eq. (1) below, which is 
simplified slightly to facilitate algebraic manipulations 
employed in developing the ensuing theory. 


The three-parameter stress-strain equations of Ram- 
berg and Osgood‘ are utilized in the following develop- 
ment with the notation altered to the o-e form. 


ASSUMPTIONS 


The analysis presented below is based upon the 
assumptions that: 

(1) Cover buckling constitutes failure. 

(2) The cover buckles as a group of parallel plates 
simply supported at the web lines and loaded edgewise. 

(3) Web buckling (in bending) occurs simultaneously 
with cover buckling. 

(4) The cover alone resists the applied bending mo- 
ment. 

(5) The secant modulus applies for plastic compres- 
sion buckling of simply supported plates. 

(6) There is no shear lag anywhere in the box beam. 


SYMBOLS 


C, = structural chord, in. (Fig. 2) 

E = Young's modulus, lbs. per sq.in. 

E, = secant modulus, Ibs. per sq.in. 

E,,, = modified secant modulus, Ibs. per sq.in 

E, = tangent modulus, lbs. per sq.in. 

K = buckling coefficient (= 3.62 for simply supported 
plates) 

L = length of box beam, in 

M = applied moment, in.Ibs. 

W = weight (compression cover plus webs) 

b, = web spacing, in. (Fig. 2) 


n = shape parameter for stress-strain curve |= 1 + 
[0.3853 /logio (o;/o2)]} [reference 4, Eq. (29)] 


p = number of bays in structural chord (Fig. 2) 
t = wing thickness, in. (Fig. 2) 
i, = compression cover thickness, in. (Fig. 2) 
te = web thickness, in. (Fig. 2) 
a = structural thickness ratio = ¢/Cs = 1/p, Eq. (1B) 
€ = axial strain, in. per in. 
= density, lbs. per cu.in. 
o, = stress on material stress-strain curve at Es = 0.7K, 


Ibs. per sq.in. 
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Fic. 2. Wing section geometry 
a, = stress on material stress-strain curve at Es = O.85E, 
Ibs. per sq.in 
Gyp = yield stress (permanent set = 0.002 in. per in.), lbs. per 
sq.in. 
Subscript 
Q = optimum design 


THEORY 


Before the optimum proportions of the multicell box 
beam can be determined, an expression must be’derived 
for the weight of the compression cover plus webs. It 
is assumed that Eq. (11) of reference 3, 


p(4p + 1) = 5(C,/t)? (1) 


may be used to determine the optimum number of webs 


required. Eq. (1) may be closely approximated by the 


relation 

ip? = 5(C,/t)? 
from which it is found that 

Pp = 1.118C,/t (la) 
or, with slight additional error, 


b= C,/t = 1/a (1b) 


Fig. 1 contains plots of Eqs. (1), (la), and (lb 
It may be seen that the weight determined from Eq 
(1b) deviates less than 1 per cent from the optimum, as 
is indicated by the dashed lines obtained from Fig. 3 of 
reference 3. Consequently, it follows that a near. 
optimum design may be obtained by dividing the wing 
cross section into square bays. 

Assuming, then, that the structural portion of the 
wing is divided into p square bays, it is a simple matter 
to determine the section weight. 

It is apparent from Fig. 2 that 


W = pL(C,t, + [1 + p]tte) 2 


From Eq. (3) of reference 3 (t./t,; = 0.4t/b,), it follows 
that ¢ = 0.4¢, for efficient design. 
Since p = C,/t, then, 


W/pL = C,t,(1.4 + 0.4a) (3 


Optimum Stress 


The three quantities required for a determination of 
the optimum bending stress are the applied and allow 
able bending stresses and the wing weight. An assump 
tion for the allowable stress, which is tenable for wings 
of transonic and supersonic aircraft, is the stipulation 
that cover buckling constitutes failure. Consequently 
the allowable stress is equal to 


o = KE,(t,/b,)? = KE,(t,/C,)*/a? (4 
using Eq. (1b), and the applied stress is 
o = M/(C&,) 
or, using a = t/C,, 
o = M/(aC,"t,) 5 


This equation follows from assumption (4). 

For optimum design, W must be a minimum. Con 
sequently, logarithmic differentiation of Eqs. (3), (4), 
and (5) will yield expressions in terms of stress and 
linear dimensions. All the dimension terms may be 
eliminated, leaving a differential equation in o and E 
which may be used to determine the optimum stress, 
since o and E, are related continuously for any struc 
tural material. 

Since the ensuing analysis applies for a constant 
value of a, the differentiation process will eliminate the 
expression in parentheses in Eq. (3). 

Logarithmic differentiation of Eqs. (3), (4), and (5 


leads to 


dW/W = (dC,/C,) + (dt,/t,) = 0 (6a) 
(da/a) + (2dC,/C,) — (dE,/E;,) — (2dt;/ts) = 0 (6b) 

(da/a) + (2dC,/C,) + (dt;/t;) = 0 (6c) 
Solution of Eqs. (6) for E, in terms of o yields the ex- 
pression 


(dE,/dc) + (3E,/c) = 0 (7) 
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OPTIMUM 


The optimum operating stress for the section may 
be found using Eq. (7) and the curve of E£, vs. o for a 
particular cover material. The graphic solution of 
Eq. (7) is depicted in Fig. 3, the use and application 
of which are similar to Fig. 1 of reference 1. 

There are two types of analytic solutions for Eq. (7). 
One involves the use of the stress-strain expressions for 
the secant and tangent moduli, and the other applies 
the three-parameter description of stress-strain curves 
in reference 4. 

Case1: Using EF, = 
the form 


a/e, E, = do/de, write Eq. (7) in 
(dE,/de)(de/do) + 3(E,/c) = 0 


However, 


or 
dE,/de = (1/e)(E, 
consequently, 


(1/e)\(E, — E,)(1/E,) + (8E,/c) = 0 


(E,/E,)(E, — E;) + 3£, = 0 
from which 


E,/E, = 4 (8) 


This expression is solved easily in graphic form by 
plotting E, and 4E, vs. o. The value of o correspond- 
ing to the intersection of these two curves is equal to 
00. 

Case 2: Using the three-parameter description of 
stress-strain curves, as is demonstrated in the Appendix 
to reference 1, Eq. (A3), 

E/E, = 1 + (3/7) (¢/o1)""! 
Thus 

E, = E/{1 + (3/7) (¢/o,)"~"] 
and, using Eq. (7) 
dE, (3/7) (n — 1) (0/0;)""! E/o 


[1 + (3/7) (6/o,)"~"]? 


do 
3E/oa 


1 + (3/7)(¢/0;)"~! 


or 
(1/7) (n — 1) (@/o1)"~" = 1 + (38/7) (¢/01)"~' 
from which it follows that 
oo/o, = [7/(n — 4) /e—) (9) 
and 
E/E» = (n — 1)/(n — 4) (9a) 
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Fic. 4. Optimum stress, modulus, and chord as functions of 


01/02. 


since a, as used in this development, actually is the opti- 


mum stress. 


Eq. (9) and (9a) are plotted in Fig. 4a in terms of 


01/02. 


A curve of 1 vs. o;/02 is included for reference. 


A comparison of the fuselage and wing optimum 
bending stresses, utilizing Eq. (A6) of reference 1 for 
the case in which the secant modulus is applicable, 


oo/o, = [7/3(n — 2)}'/"- 
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Nonoptimum design relationships. 


Fic. 5. 


reveals that the ratio between the wing optimum stress, 
oo, and the fuselage optimum stress, oso, is 
‘ ‘ 1/( l 
Fw/ Fn = [3B(n — 2)/(n — 4)" ” (10) 
For most structural materials, and aluminum alloys 
in particular, nm > 6; consequently, o,0/o¢ generally will 
exceed unity by only a few per cent. 


Optimum Proportions 


When oo has been determined, the optimum struc- 
tural chord may be found from Eqs. (+) and (5): 
From Eq. (4), 


t, = aC,(c/KE,)” 
and from Eq. (5), 
i= iM ‘(aaC,”) 
By eliminating /,, 
C3 = (M/a?)(KE,/o*)"" (11) 
and, therefore, 
t,®> = Ma/KE, (12) 
Upon substitution of C, and ¢, in Eq. (3), an expres- 
sion is obtained for the weight of the compression cover 
plus webs: 
W = pL(1.4 + 04a) (M4/[a?2?KE,o* |)" (13) 
The values of Cyo, to, and Wo are obtained by using 
oy and Fin Eqs. (11), (12), and (13) instead of o and 
E,. 
The ratio of a nonoptimum chord size to an optimum 


chord size for any J/ and @ may be found from Eq. 
(11): 
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&, Cs = (FE,a03 Eo)’ . (lta 


Similarly, for optimum and nonoptimum weights, 


using Eq. (13), 
W/Wo = (Eov®/E,o*) l4b 


Multiplying Eq. (l4a) by Eq. (14b) leads to the re 
sult, 


(C,/Cw) (W/W) (¢/o0) = 1 15 


The identical expression involving shell radii instead 
of beam chords has been found to exist for long circular 
cylinders in pure bending [reference 1, Eq. (11)]. 

The variations of ¢/oo and IV/W» caused by depar 
tures of C,/Cyo from unity are depicted in Fig. 5, which 
was constructed in a fashion similar to that used for 


Fig. 5 of reference 1. 


Simplified Design 

A close approximation may be obtained to the opti 
mum structural chord through the use of a simplified 
design formula derived as follows: Rewrite Eq. (11 
to read 


Ce os (M?K Eo atgy*)' . 


Since 


then, 


ee" ; ” | K E (“) | sc Tr 
(M?E/a‘to,*) ’° Eyo\oy 


The right side of Eq. (16) is a monotonic decreasing 
function of o;/o2. When Eq. (16) is plotted as shown 
in Fig. 4b, a value for this function may be found which 
vields a minimum weight deviation for the indicated 
range of o;/a2. Since the cover is assumed to be simply 
supported, A = 3.62. By trial and error, using Fig. 5, 
the value of the left-hand side of Eq. (16) corresponding 
to a minimum weight deviation is found to be 1.167, or 
7/6. Thus, a close design may be obtained by use of 


the expression 
bi - ’ om a\ 1/6 joi 
Cxo,) = (7/6) (.M*E/a4o,') (17 


The curve for Cy.,)/Cyo is plotted on Fig. 4b, from 
which it may be seen (using Fig. 5) that the maximum 
weight deviation is less than 1 per cent for 8 << n < 20 
For n = 6 the error is about 3!/2 per cent. However, 
the value of for most of the common aircraft struc 
tural alloys lies between 8 and 20, and, consequently, 
use of Eq. (17) for design will lead to a negligible weight 
increase. 

Since the value of o,, usually is known for a particular 
alloy and since the stress-strain curve must be avail- 
able to find o;, it would be desirable to use the yield 
stress in Eq. (17) instead of o,. This may be done, 
since the additional error is not serious, as will be 
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Table 1. Chord, Weight and Stress Errors, for Various Structural Materiais, Incurred 
by Use of the Simplified Design Procedure 
Aluminum Alloy Sheet Stainless Steel Sheet 
17S-T 24S-T* 75S-T* Stainless We Type 301¢ Type 304* 
[tte Room Temp. | Room Temp. 400°F | Room Temp. 400°F Room Temp. Room Temp. | Room Temp. 
[ € (108 psi) 10.67 10.50 9.60 10.42 9.35 28. 30 26.00 25.68 
| a, (107 psi) 40.10 45.60 35. 30 72.80 31.00 191.0 33.10 30.60 
| a, (107 psi) 34.80 41.25 32.30 69.10 28.25 171.0 30.10 25. 20 
| %» (107 psi) 41.40 46.00 36. 00 71.50 31.75 185.0 37.05 36.70 
| ola, 1.152 1. 106 1.093 1.053 1. 088 1.117 1.100 1,213 
| ao 7.22 9.83 10.95 18.10 11.50 9.08 10. 30 §.55 
Telo, 4. 133 1.021 1. 001 - 960 -994 1.040 1.011 1.372 
o,/0,5 . 968 991 - 980 1.018 975 1.032 - 893 834 
Cio, )/ aio) - 984 - 996 - 990 1.009 - 987 1.016 - 945 -913 
| Cy (2 ,)/Cso 1.113 . Be - 997 -951 -990 1.037 1. 008 - 880 
Cy (o,1/Cno 1.096 1.013 - 986 . 959 977 1.052 - 954 . 804 
(W/W,)-1, % 0.9 0.1 <0.1 0.4 <0.1 0.4 0.3 1.4 
| oo, .875 . 985 1.015 1.040 1.025 . 940 1.045 1.125 
| ole, ,-1, % -3.9 -0.3 -0.6 +1.8 -0.6 +0.9 -5.7 +28.8 











*see Table 1 of Reference 1 for sources of material data 


lemonstrated: 
Csayp)/ Cs0 a [Csceyp) Cs¢o1)] [Cs Co] 
Values for C5¢4,)/Cyo may be found from Fig. 4b, while 
Otis C, a ‘iti (o, ‘Oyp) ; 


from Eq. (11). Consequently, Cy¢,,)/Cso may be com- 
puted for any material, and the corresponding values of 
W/W) may be found. This has been done in Table 1, 
from which it is seen that the weight deviation has a 
maximum value of 1.4 per cent, while most of the errors 
ire less than '/, per cent or are negligible. Therefore 
itis feasible to use the following design expression : 
7/6) [APE/ate,,*}' 


Coo = (1S) 


rhe relation between the actual section stress and 
the optimum stress (¢/o,9), as well as the discrepancy 
between o and gy», is included in Table 1. The plus 
ind minus signs indicate section stresses respectively 
The most serious 
This 


is expected because the constant in Eq. (18) was se- 


greater and less than the yield stress. 
discrepancies occur for type 304 stainless steel. 


lected for materials with m lying between 8 and 20. 
DISCUSSION 


Choice of Modulus 


The theoretical investigation of plastic buckling of 
plates reported by Stowell? provides expressions for 


the various effective moduli pertinent toedgewise loaded 
plates with various types of edge restraint. The modu- 
lus obtained for simply supported plates, 


E,, = E,{(1/2) + (1/4)V1 + (3E,/E,)] (19) 


“3m 


is slightly lower than the secant modulus, obtained 
for plates simply supported along one unloaded edge 
with the opposite edge free. Since the cover of a 
multicell box will probably behave like a series of simply 
supported plates arranged side by side with the buckle 
ridges and valleys alternating from bay to bay, it would 
be of interest to estimate the difference in weight ob- 
secant (E£,_) 


m 


tained from the secant and modified 
theories. 

From Eq. (8), /i/E, = 
ating at the optimum stress according to the secant 
Substituting this value into Eq. (19), 


This assumes that the 


'/, for all materials oper- 


modulus theory. 
E,,,/E; is found to be 0.831. 
optimum stress according to the /,,, theory does not 
depart appreciably from the F, theory. 

Since the weight is proportional to the product of the 
structural chord and cover thickness [Eq. (3)], then 


from Eqs. (11) and (12), 
W~ Ct, ~ E, 


or 


It was assumed implicitly, in the analysis above, that 
the secant modulus is applicable to plastic bending of 
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the webs. It must be assumed, for purposes of this 
comparison, that £,,, applies to cover and web alike. 
This is consistent although not necessarily correct. 

It follows that the ratio of the E,,, weight to the F, 
weight equals 


(FE, /E,)~"* = (0.831)7'* = (1.032) 


Sm 


Thus, the secant modulus theory will yield a design 3.2 
per cent lighter than that employing the modified secant 
This problem (of the correct modulus) war- 
both through tests and 


modulus. 
rants further investigation 
more extensive analysis. 


Rectangular Bays 

The deviation of the bay from a square will not 
affect the expression for optimum stress, provided the 
bay proportions are maintained constant for all values 
of C, and ¢t, although the computed weight will change. 
This change will be small, however, as is demonstrated 
in Fig. 1, if the change in shape is not severe. This 
follows from the solidity curves of reference 3, which 
demonstrate the small change in weight which accom- 
panies a change in the number of bays for a structure 
designed according to Eq. (1) above, and from the con- 
sideration that a deviation in bay shape from a square 
to a rectangle may be considered as an increase or a 
reduction in the number of bays. The actual magni- 
tude of this effect must be determined separately for 
each particular case, especially when the bay aspect 
ratio departs radically from 1. 


Shear Effects 


The previous analysis must be modified to in- 
clude the effects of shear before it can be used for gen- 
eral wing design. Although this modification might 
not be severe for some types of loading and for a few of 
the many possible wing shapes, nevertheless, the general 
case must be investigated before any extensive conclu- 


sions can be drawn. 


Shape Effects 


The variety of wing plan forms and methods of root 
attachment being used, and investigated for use, in 
high-speed aircraft requires examination for their in- 
fluence upon the results of this investigation. The 
effects of sweep and taper alone may require large 
structural weight increases. Furthermore, the applica- 
tion may demand heavier wings from. a performance 
standpoint alone. For example, a guided missile 
might be required to maneuver up or down with equal 
facility, in which case both covers would be designed 
identically for compression loads. Under these cir- 
cumstances the expression for the total structural wing 


weight would become 
W = pL(2.4 + 0.4a)C,t, 


Since missile wings are extremely thin, the second 


term in the parentheses is negligible compared to the 
first, and the structural wing weight becomes 


W = 2.4pLC,t, 


If the root connection does not develop the full cover 
load because of riveting or fitting difficulties, it will be 
necessary to increase the cover and web gages for some 
distance outboard of the root in order to transfer the 
shear and moment into the fuselage. The obvious re 


sult is an increase in weight. 


Simplified Design 

An examination of Table 1 reveals that appreciabk 
deviations from the yield stress result in small weight 
errors based upon optimum design. The largest errors, 
for type 304 stainless steel, were 1.4 per cent in weight 
and 28.8 per cent in stress. 

Assuming a 10 per cent weight error for a material 
with m = 10, which is a typical value for the materials 
listed in Table 1, it is found (using Fig. 5) that C,/Cy = 
1.245 and o/o,, = 0.71. In other words, 29 per cent 
deviation from the yield stress results in a 10 per cent 
weight error. 

The results of the analysis leading to the simplified 
design procedure infer that a box beam designed to 
operate in the neighborhood of the yield stress will 
weigh close to the minimum value as is indicated by 
Eqs. (17) and (18), which show the optimum chord to 
be proportional to o “?, and by Fig. 5, which shows 
small weight errors attending appreciable discrepancies 
in C, when compared with Cy. If shear effects and 
shape factors are small in a multicell wing designed for 
buckling near yield, it will be designed close to opti 
mum, and little weight reduction will be obtained by 
attempts to improve the design in order to approach the 
yield stress more closely. This means, for instance 
that it may not be necessary to use tapered sheet in some 
cases that apparently require it, since the standard 
gage variations may permit a type of construction ata 
weight nearly equal to that employing the tapered sheet 

In the preceding analysis and discussion the emphasis 
has been upon percentage weight change instead ol 
actual poundage. Obviously, in a heavy wing a few 
per cent will equal many pounds. Thus, the utility ol 
the included information is a function of a particular 
case to which it may be applied. 
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The Stress Analysis of Arbitrarily Shaped 
Wing Structures with Shear Stressed 
Skin and Walls 


H. F. MICHIELSEN* 
Fokker Aircraft Factories 


SUMMARY 


The internal load of a statically indeterminate wing structure 
an be decomposed into a basic system and a number of self- 
equilibrating load systems. Each of the latter involves two or 
three rib segments and three or two spar segments, as well as the 
two upper and two lower panels bounded by them. The calcula- 
tion of the forces in these statically determinate systems is pre- 
sented. The fundamental problem is the calculation of the mul- 
tiplying factors (statically indeterminate quantities) of the 
self-equilibrating systems. 

With the aid of Castigliano’s Theorem of Least Work, the 
redundant quantities are found as the roots of a system of linear 
A suitable method of solution, first suggested by 
Gauss, is presented. The coefficients in these systems of linear 
equations are obtained as the integrals of the products of two self- 
equilibrating systems or of one such a system and the basic load 
The loads can be decomposed into one to three unit 


equations 


system. 
loads, depending upon the structural element, and the integra- 
tion need be carried out only for these unit loads. 

The coefficients of the equations are then obtained through a 
summation of the integrals of the products of the unit loads, after 
each integral is multiplied by a suitable factor. The main prob- 
lem is, therefore, the systematic determination of these multi- 
plying factors. Details are given in a number of tables. 

A special section deals with the evaluation of the integrals of 
the products of the unit loads. However, this part of the calcu- 
lations is somewhat uncertain because of the unreliability of some 
of the assumptions. 

A concluding remark indicates the possibility of an extension of 
the procedure into the plastic region. 


(1) INTRODUCTION 


pe METHODS are known by which the stress 
distribution in wing structures can be calculated 
more or less adequately. Wing structures are gener- 
ally statically indeterminate to a high degree and the 
essential problem is mainly to determine how much ef- 
fect the torsional rigidity of some elements and the 
bending rigidity of some combinations of elements has 
im resisting the applied torque. Earlier investigators 
strived at finding analytical solutions in a closed form. 
They had to assume continuous and analytically defined 
variations in the dimensions and rigidities. Among 
other things the transverse stiffening elements (ribs) 
had to be distributed continuously. 

The development of modern multiengined airplanes 
with cantilever wings, containing gas tanks, retractable 
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landing gears, and occasionally luggage compartments, 
made it impossible to maintain the assumption of con- 
tinuously distributed transverse stiffeners. As a con- 
sequence, the solution in a closed analytical form be- 
came useless. It was replaced by calculations based 
on a finite number of redundant internal forces or 
moments. The actual number of ribs, defining the 
degree of redundancy, is far too large for convenient 
calculation. On the other hand, the replacement of 
a number of ribs by a single substitute transverse stif- 
fener is usually an assumption accurate enough for prac- 
tical purposes, so far as the calculation of the redundant 
quantities is concerned. After a solution is obtained, the 
loads, calculated for the substitute rib, can easily be 
distributed among the actual ribs in accordance with 
their individual stiffnesses. 

The fundamental limitations of the method described 
in this paper are: 

(1) The wing is composed of an arbitrary number of 
spars, transverse stiffeners (ribs), and an upper and a 
lower set of panels of cover sheet. 

(2) None of these structural elements is able to carry 
torsion by itself. 

(3) The cover sheet is subjected only to shear stresses 
and possibly stresses arising from a distributed trans- 
verse load. Its edges are not loaded perpendicularly 
to the edges in the plane of the cover sheet. 

In the calculations of Section 4 it appears advanta- 
geous to make a few more assumptions. They are made 
in order to simplify the calculations of the statically 
determinate simple cells. The most general solution of 
the simple cell problem constitutes some further dif- 
ficulties not essential to the procedure outlined in this 


paper. 


SYMBOLS 
A = area of cells 
Aa = area of panel between connecting rivets 
Ay = area of flanges 
An = mean area of cells 
Aw = area of projection of spar web on Y-Z-plane 
Awa = real area of spar webs between connecting 
rivets 
E = Young’s modulus 
F,, Fs, Fe = unit loads on flanges of spar webs 
Fac, Fe = SF.dS, J FFAS, ete. 
G = modulus of rigidity 
I = moment of inertia 
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= load 
Lo = basic load 
L;, L; = self-equilibrating load systems 
M = bending moment 
MM, = torsional moment 
N = normal load in spar flange 
P, = ,unit load on pair of panels 
Pram = Sf PpdS 
R.. Ry X = unit loads on rib segment 
Raa, Ri = SRdS, J R:RAS, ete 
S = shear force 
dS = dx/EI, dx/GA, dxdy/Gt 
Wa, Ws, W. = Unit loads on web of spar segment 
Was, Wre = fWidsS, JWiWAS, ete. 


X = redundant quantity 
= spar distance in X-direction 


Aa = width of projection of spar segment on X-Z 
plane 

b = rib distance in }-direction 

f = spanwise numbering ( }-direction) 

g = chordwise numbering (X-direction) 

h = mean height of rib segment 

he; = effective height of rib segment at station p 

hap = height of rib segment between connecting 
rivets at station p 

hp = total height of rib segment at station p 

4,7 = numbering of self-equilibrating load systems 

l = length of spar segment 

m = numbering of sections of spar and rib segments 

p = numbering of stations of rib segments 

n = number of sections of spar segment 

n = degree of redundancy 

q shear flow 

qm = mean shear flow of a cell 

dp = shear flow in a panel 

Gr = shear flow along the circumference of a rib 

t = thickness of panel or web 

a = degree of taper in a 

8 = degree of taper in h 

by = he/h ‘ ( spar segment 

6 = degree of taper in h, 

n = degree of decrease in E 

rN = nondimensional spanwise coordinate 

I = front spar of a cell 

II = rear spar of a cell 


(2) Basic CONSIDERATIONS 


The fundamental problem is the determination of 
the statically indeterminate quantities and, in sequence, 
that of the internal moments, shear forces, and shear 
flows, to be referred to as internal loads. The calcula- 
tion of the stresses will be discussed only as far as it is of 
importance in the determination of the internal loads 
(Section 9). The are 
assumed to be known. 

The rigidities of the individual structural elements 
are assumed to be vanishingly small for types of loading 
to which they are not subjected to an appreciable de- 
For this reason spars and ribs are assumed to be 


external loads and reactions 


gree. 
infinitely weak in torsion, and cover sheets are assumed 
to be subjected to shear alone. The assumptions made 
result in from six to about 20 redundant quantities in 
most airplane wings. Although this number is large, 
the computation involved does not present undue dif- 
ficulties, as will be discussed in Section 3. 
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In a statically indeterminate structure, the interna] 
loads L, acting upon the structural elements. can be 
represented by the expression 
Xe (] 
where Ly is an arbitrary system of internal loads that 
balances the known external loads, L; is a system of self 
equilibrating internal loads, and n is the number oj 
independent self-equilibrating systems L; which can 
exist in the structure. Finally, Y; is a statically in 
determinate factor. Eq. (1), with arbitrary factors 
\,, represents all the possible internal load systenis that 
can balance the known external The 
values of the factors XY; must then be calculated fron 


loads. actual 
the requirements of consistent deformations. 

It should be noted that, in the calculation of the 
systems Ly and L,, the sectional properties of the struc- 
tural elements are not needed. The geometry of the 
structure suffices, together with the magnitudes of the 
external loads as far as Lo is concerned. On the other 
hand, the values of the X; can be calculated only if the 
sectional properties are known. 

The calculation of the loads Ly and L; can be carried 
out with any degree of accuracy, irrespective of what 
the material properties and the numerical values of the 
actual loading are. In contrast, the factors X, depend 
greatly upon the moduli of elasticity, which vary con 
siderably with the loads and which are somewhat un 
certain. It therefore seems of no use to strive after an 
accuracy greater than, for instance, +5 per cent in 
evaluating the definite integrals /L,L,dS, to be dis 
cussed in Section 3 and used in the computation of Y 

Finally, it should be noted that the load systems L 
and Z; are not uniquely determined. This can be seen 
from the following development: 


L = Lo + XiLy + Xele + X3h3 = (Lo + Vili + 
Y3L3) + (X; — Vi;)Li + Xole + (X3 — Y3)L; = 
(Lo + Vili + Yils) + (X%1 — ¥1 + X2)Li + X2(L2 - 


ibe) : (Xe = Y3)L3 = 8 * a b fig & = Be Fg oe Xe 


Ihe load systems Lo*, L;, L.*, and LZ; can be obtained 
from the load systems Lo, L;, L2, and L; by means of 
linear transformations. It is easy to prove that all the 
requirements mentioned before are satisfied by the new 
systems as well. The most suitable choice of load 
systems will be discussed later (Sections 3 and 7 


(3) FUNDAMENTAL SOLUTION 


Before the loads Ly and LZ, are calculated, the general 
solution of the statically indeterminate problem will be 
Use will be made of Castigliano’s Theorem 
The strain energy will be calculated 


discussed. 
of Least Work. 
by integration of the product of load and distortion 
each structural element. As is well known, the strait 
energy dIV can be expressed as: 


dW = (1/2)L*dS (2 








pos 


abl 


internal 
can be 


ads that 
n of self 
mber of 
iich can 
cally in 

factors 
21s that 
e actual 


ed fron 


1 of the 
le struc- 
y of the 
Ss of the 
he other 
ly if the 


> carried 
of what 
2s of the 
depend 
ary con 
rhat un 
after an 
cent in 
be dis- 
n of X, 
stems L 
be seen 


Y,L, + 
V3) L3 = 
Yo(Le - 
+ X;*L, 
»btained 
neans Ol 
t all the 
the new 
of load 
L 7 


- general 
1 will be 
[Theorem 
Iculated 
yrtion 10 
ie straif 








STRESS ANALYSIS 
where L is some internal load and dS the stiffness factor 
of the structural element. 

It may be noted that the dimension of L varies. 
instance, L can be a moment, a force, or a shear flow. 
The corresponding dimensions of dS differ also, for 


For 


instance, 
dx/EI = (force! X length~') 
dx/GA (force~! X length) 
dx dy/Gt = (force! X length*) 
The dimension of IV” is, of course, always (force X 
length). 
Eq. (2) is correct only if any deformation is caused 


exclusively by the internal load that appears in the 
equation involved. This is true in the case of spars of 
variable height, provided the shear force S is meant to 
be the part of that force acting in the web. In general, 
the wing structure must therefore be imagined as con- 
sisting of bars subjected to axial force, such as the 
flanges of spars and ribs and possibly vertical stiffeners, 
and plate elements subjected to shear alone, such as the 
webs of spars and ribs and the sheet covering. 

After substitution of Eq. (1), integration of Eq. (2) 
leads to 

W= (1 2» f (lo + >> XiL,)%dS 
i=1 

The partial derivatives of IV with respect to any of 

the factors Y; must vanish: 


OW, ox; = 


/ (lo+ > XL)LAS = 0 
i=] 


or in another form: 


[ tetas +> xf L,LAaS = 0 (3) 
. i=1 


Since 7 can assume any value from 1 to 7, » equations 
of the form of Eq. (3) are obtained with » unknowns 
X, to Xp. 

The present system of simultaneous liner equations 
has certain properties that make its solution practicable. 
These properties are: 

(1) The coefficient of X; in the jth equation 
(f'L,L,dS) is identical with the coefficient of X; in the 
ith equation (f'L,L.dS). 

(2) The coefficient of X; in the jth equation (f'L,L; 
dS) is generally much smaller than the geometrical 
mean of the coefficient of X; in the ith equation 
(SL dS) and the coefficient of XY, in the jth equation 
(SL dS). 

In order to take advantage of the latter property, it 
is advisable to choose the internal load systems L; so 
as to have the least amount of overlapping. Conse- 
quently, each load system should influence the smallest 
possible number of structural elements. 

Systems of linear equations of this kind are prefer- 
ably solved by a method suggested by Gauss.'! Since 
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Fic. 1. Plan view of part of wing structure. 


the method is not well known, it is shown in the Ap 
pendix in a form modified slightly by the author 


(4) THE SELF-EQUILIBRATING SYSTEMS 


Fig. 1 represents a portion of a wing structure. Three 
spars divide the wing into four chordwise sections, each 
of which is subdivided by ribs in cells. A cell is shown 
in the picture by shading. 

Two adjacent cells have either a spar segment or a 
rib segment in common. Under the assumptions made 
in regard to rigidities, the common faces of the two 
cells are coincident if their four corners are coincident. 
The forces necessary for deforming the two cells until 
the four pairs of corners are in coincidence give rise to 
stresses in both adjacent cells exclusively. Except for 
a multiplying factor, these stresses are determined 
unambiguously by the geometry of the cells. 

The number of redundant quantities is, therefore, 
basically equal to the number of common faces of the 
cells, although the faces marked 0 in Fig. 1 are not to 
be counted so long as there are no spars in the leading 
edge and the trailing edge. 

As yet no assumptions have been made regarding the 
shape of the individual cells. A few will now be made 
in order to reduce the work involved in calculating the 
forces in the cells. In most cases these assumptions 
are not in contradiction with common structural usage. 
If they are, the internal forces must be calculated in a 
different manner, but the general procedure of this 
paper can still be followed. The assumptions are: 

(1) All spar and rib segments are plane, or, in other 
words, their two flanges and web are contained in one 
plane. In addition, this plane is assumed to be parallel 
to axis Zin Fig. 2. The spars may have a break at each 
rib, while the ribs are straight and parallel to the X-Z- 
plane. 
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\+Z 
qd, 
Fic. 2. Cell under torsion 
(2) The panels of sheet covering are ruled surfaces, 


the generators passing through the intersection of the 
adjacent spar planes. The intersections of these sur- 
faces with the planes of the spars are therefore straight 
lines, in some cases skew, while the intersections with 
the planes of the ribs are arbitrary. 

(3) The flanges of the spars are straight between 
adjacent ribs but they may be broken at the ribs. 

On the basis of these assumptions the following ex- 
pressions can be derived for the internal loads. The 
derivations of these expressions are omitted in order to 
save space. 


(a J oe in a Cell Bounded by Two Parallel Ribs (See 
ig. 2) 


If the average rib heights are denoted by the symbols 
h, and hs in accordance with the equations 


hy = A;/a, hs = A;/ds (4) 
then the average shear flow is 
dm = Me (ayhz — Azh) (5) 


With the aid of gm2 the local shear flows in the panels 
can be calculated: 


(6) 


dp = Qm2(a1d3/a") 
In particular, 
(7) 
It is now convenient to introduce the nondimensional 
quantities a2, 82, and rx». a», and 82 represent the taper 
in width and height of the cell, respectively, while ). 
is a coordinate parallel to the Y-axis having the limits 
—land +1. 


Upt = Qm2(Gs/a1), Jp3 = Ym2(a1/ds) 


@ = a2(1 — reas), G1 = G2(1 + a2), a3 = a2(1 = 
h = ho(1 — roBe), fr = ho(1 + Be), hs = he(1 — Bo) 
(8) 
With 
2A2 = ayhz + a3hy = 2a2h2(1 - a2Be2) (9) 
yields 
dma = Mp/2A; (10) 
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and 


1 — & 
dpi = Im2 ii“ 


5 tl (11 
dps = m2 oer 
—_ 


The bending moments in the spars can be calculated 
from these shear flows. One finds, for the XY compo. 
nents of the bending moments in spars I and ITI (see 


— Bef 1—*\_ 

20 (= ~) 

witn(*)f p22 # AE) a 
a2/L4(1 — aeBe)_J\1 — Arcae 


The equation shows that the moments in the two 
spars are equal and opposite, since the external moment 
M, is equal to zero for the entire cell (pure torsion), 
Moreover, the bending moments vanish in the spars of 
conical cells or rather in cells whose width and average 
height are proportional (a: = £2), although the spars 
may have a different taper. Finally, the bending 
moments vanish at the end ribs (Ag = +1), while their 
variation is substantially parabolic (1 — ,”). The 





Fig. 2): 


Maru = +Qmabohe = 





multiplying factor (1 — A,a2) in the denominator causes | 
a slight deviation from the parabolic form. 

Since the flanges of the spars were assumed to be 
straight, the effective spar height 4,—that is, the dis- } 





tance between the centroids of the flanges measured in 
the Z direction—can be expressed in the following man- 
ner: 
i har = he(yer + 5er ) 
hor = he(yor — Aeder ) Ju 72 
Vist = ha(yer — der) 13 
(15 
jj = h(¥: dor1) , 
Matt to(yort + Sort 
han = he(yer1 — Arde) Me 
Ua = Neyer — don 
The Y components of the forces in the upper flanges 


follow from Eq. (12): 
Nain = 
1 _— do? 


F mabe (* = *\) | (14 
2 (L — Asa@e)(Yer,11 — Acder,11) 


This last equation permits the determination of the 
shear flows in the spar web: 


a 


j 1 -_ Qty” _ 
((1 — Araz)? 
l- a2? 


a2 = Bo | ae. 
ayyor,1t — Se1,11 L(1 — Azae)? 


= + m2 . 


Yorm1” — 921,11" } 
(Yer,1r — Agder,11)? 

(15 
It is noted that Eq. (15) is considerably simplified when 
v2 = 1 and 6. = fs, which is true, for instance, for cells 
with rectangular cross section and flanges in the cor 
ners: 


Quit = * Gm [(1 — B2*)/(1 — A2B2)?] 
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This equation is identical with Eq. (11) if a, is replaced 
by 8 in the latter. 


The loading of the ribs depends upon the manner in 
which the torque 1/, is applied. Two cases can be dis- 
tinguished. 

(b) Two Adjacent Cells Have a Common Rib Segment (See 

Fig. 3) 

The self-equilibrating system is defined in such a 
manner that in the rear spar (II) the bending moment 
at the middle rib (3) has an X¥ component 


Maun = +1 (16) 
In the front spar (1) one finds 
Maar = —1 (17) 
and, since the bending moments decrease to zero at the 
end ribs (ribs 1 and 5) , 
Mastin = () (18) 
one obtains 
Sa = =1/b, Sian = 1/4 (19) 
fhe corresponding reactions at ribs 1 and 5 are trans- 
mitted to cells 2 and 4 in the form of torques. 
Me = +4;/b,, Mu = — as/b5 (20) 
When there is a break in one or two spars at the mid- 


dle rib, the component of the bending moment in the 
spars causes bending in the middle rib: 


Ader Adgi11 ‘ 
Mat. = =( be 4 b, ) (21) 


c) Two Adjacent Cells Have a Common Spar Segment 
(See Fig. 4) 
The self-equilibrating system is specified in such a 
manner that 


Mins = +] (22) 
From this follows 
Me = —Mp* = —- 1 (23) 


These equations suffice to calculate the loading of the 
individual structural elements. 


(5) THe Unit Loaps 


The definite integrals {"L,L,dS, derived in Eq. (3), 
are conveniently calculated with the aid of a small num- 
ber of unit loads which can be combined into the actual 
load systems L;. The integration has then to be carried 
out only for these unit loads. The values obtained 
must be multiplied by factors that depend only upon 
the geometry of the structure, and the necessary com- 
putations can be easily made in a tabular form. 
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Fic. 3. Two cells having common rib segment. 
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Fic. 4. Two cells having common spar segment 
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Fic. 5. Unit load on upper panel. 
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Fic, 6. Unit loads on rib segment 
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Fic. 7. Unit bending loads on spar segment 


(a) Concerning the Panels (See Fig. 5) 


Since it was assumed that shear stresses alone can act 
along the edges of the panels, no other load systems 
need be taken into account. The unit load is defined 
as that corresponding to a unit torque eventually acting 
upon the cell whose structural element the panel is. 
Consequently, the unit load, to be denoted by the 
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symbol P,,, causes a shear flow according to Eq. (10): 
Gm = +1/(2An) (24 


The direction of the shear flow in the upper panel js 
shown in Fig. 5, and the value of A,, is given in Eq. (9), 

For completeness it is noted that small edge loads, 
perpendicular to the panel, are necessary for equilibrium 
when the edges that are attached to the spars are skew 
These perpendicular loads give rise to small bending 
moments in the panels, whose contribution to the strain 
energy may be negligible. These bending moments are 
proportional to the shear flow and thus can be included 


into the unit load. 


b) Concerning the Ribs (See Fig. 6 

The rib loads, derived in Section 4, consist of a con- 
stant shear flow along the rib flanges and of shear forces 
and moments in the end sections of the rib segment, 
The distribution of the shear force along the height of 
the rib influences only the vertical stiffeners. As the 
strain energy stored in these vertical stiffeners is small 
and therefore neglected, the loads acting upon the rib 
can always be decomposed into the three unit loads 
shown in Fig. 6. Each of the three unit loads R,, R,, 
and R, corresponds to a negative unit moment, applied 
to the rib in three different ways, and causes a shear 
flow along the entire circumference of the rib segment 
amounting to 

gz = +1/(24,) (25 

A, must be calculated according to Eq. (4). 


(c) Concerning the Spar Flanges (See Fig. 7) 


The bending of the spars can also be decomposed into 
three unit loads, two of which, F, and F,, show a linear 
distribution of 1/7, from 1 to 0 and from 0 to 1, respec- 
tively, and the third, F;, a parabolic distribution of M, 
from 0 through 1 to 0. It is to be noted that Eq. (12 
represents a parabolic distribution in a first approxima- 
tion only. In Fig. 7 the actual moment distribution 1s 
shown dotted. 
means of the parabola has the advantage of reducing the 


number of the unit loads from four to three and seems 


The approximate representation by 


permissible, since the deviations change sign and since 
the multiplying factor according to Eqs. (12), (20), 
and (23) is small. It is of the order of magnitude ol 


(a — B)/2. 


d) Concerning the Spar Webs (See Fig. 8 


The loading of the spar webs can also be decomposed 
into three unit loads. The loading W», corresponds to 
the loadings F, and F, shown in Fig. 7; yet the mag- 
nitude of the unit load is different. From 


mo = 1/(2Aw) = 1/[b(ha + Aes) (26 


where A,, is the projection of the web surface upon the 
YZ-plane, measured between the centroid lines of the 
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J. (10): flanges, there follows, because of Eq. (7), 
(24 761 Tel 
‘ d3n = YUmp = f f 
h. a a) IP —— a oe £¢g. 
panel is ” ma Via + he L, 77 | Gm fPAlap (6)p> ¥ 
Eq. (9). On the other hand, load system F, with 1/,, = +1 a 
e loads, and \/,, = 0 corresponds to a shear flow 32 
librium 1 / (Bh 12 bp (10)p5 (20) 
Oa)! = (Dil) | 
re skew . ” | 
: a | 
bending Hence, ef f | 3 
le strain : } (2 J 
h. ceae hes . 
ents are (an) WV, = l (qm) Fe, = ( Jaa )F. 23 +f 
ncluded ha + hes ha + he a 
7 
(27) 32| + a | *(G)5| (20) 
[he unit loads IV, and IV, correspond to the unit 2 
load F, of Fig. 7, which latter really represents two 
f a con- cases of loading. Here again the magnitude of the 
ir forces unit load is different, since the moments J/»2* = +1 7a@4le 2 Fuib segment 32 
egment. and WJ. = —1, by which the loads W, and IW, are de- 2. ali. fir nee ae f Alg=-/| Fo. 
; : ; : 7 Ja %= pw 77 yy ree  ( 3 
eight of fined, yield a bending moment in the center section of o| 2A,2| © 2As2 2-M)32-(e)32| “\"" 2Abp 
As the the spar segment, determined by Eq. (12), instead of | 42 i + fue ce + 5 )po) - i) (20) 
is small unity as assumed in the unit load F,. Ass Ye 
the rib It is observed that the loads WV’, and HW’, are about af +f * Aes Q32 -1\+/7bo (7) /e3)\ 
it loads conform to the load IV’,, as will be shown in Section S 23 - Fe ee +f '-/Tbo =f (10) { 
R,, R, [Eq. (46)]. Therefore, the loads WW, and W, can be | pd | 
. = ee Sa ae a ae *) 4 472% %e 
applied replaced approximate ly by a »(A u A =) end a » x *) ’ e ae Fe De " (7) 
a shear A,/A»), respectively. Advantage of this statement | 32 Pe Pill -lt7 132 -) (10) (20) 
~ i aoe ; » 7 ~ 52! 0 2 oan ¥ 
segment was taken in Table 10. | Br) -(22, i220 Z, - (24) +(23), (20) 
| Aan 0. 
_ =a ae fi \- ) — 
(95) — al De Ap cp ”* ("he ad (23) 
a Aap 
- |Z < —we 5? _f |-/7) 
W, fr ¥ +f “ "Aue [a2 The ~ 
TOM I, = f 9a Isp ._ V7) po 
2 [iy~4 se| . oy Mla M10) ©! 
2 -_ V4 ane or both of the spars has a break of points 3.4 
a linear 
r™ andfjor 3.3, the factor far R, becomes:-[4 et _ 9 Gt 
— 2 
respec- 4p &, 
n of M, and the factor for A. becomes: + 4023 _ ae) 
3 9) 2 eg 
aq. (12 aceoraing to £y. (Pt). The factor for Fe changes in such 
i — amanner, that the factors for AL, Fy and Re ada 
ution 18 — 4° to the value, given {Gr ki in the table. | 
f10n by 
cing the 
d seems Table 3 Flanges in spar segment 23 
- ~ —— 7 
” £: [Fale] 5 | Mea= et AZZ Ee 
ae (2 }y ] — — 4 
asic 22 _. %22-/Ip ; | 
itude ol af | - | a9eeig Bee Ae + (32po - (12), (23) 
/- %2)(822 -/322) 
12 +* * a 36, _ ] 
| | W220 Aee) Flee 6) 20) 
|_ 4+ A%z2NOp2-Az2) _ " (2) 
” -_ 4('- ee ee) 5)>2 f 
imposed ~-2%e_,_AeSi22 | _ p39) 
+ lee 2 
onds to 23) = % ‘Sel %2/%22) _ (2) ,/23) 
a __ 242 %4-/924 -(38) ' 
he mag fz Z 0+ Cay Yll-O24 fey) °4 | 
(1- %24)Op4 -/324) 
14| - . - (3%, = 
“s Ce C7. f 41-0124 (24) vee (8) ro | 
(20) (1+ O24) (Az4-A24) |, (35, fey" 
34 ie + +(35) ‘a 
the £9. (15) 4(1-A2, (324) a | 
1 2 ae 
por ; ; 25 a +s “e. rar "~ + (92) = wanes) 
s of the Fic. 8. Unit shear loads on spar segment. il e6/"2: 
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‘Mt ic. Table S Series of cells - fo] 
tela (4) = ” 6 
f rth f 3 s 
| ie ic, 
> /) a a a3 Qs } 
| /2) A A, | As an 
L213—1.4 (3) | Mpg “(2lros ahs | ahs | 
* . be) | Mpls ss | = 
(5) 4) -/4) G3 - Ay | %hG- Ash; 
ei @2) 245 16)| 3) +/+) \oh 4a;hy| Ghgrash; 
= , 3 2 TA 2. i 
qr] ae [ae taeetee| | 
~ (3) ye) 4p b, 
Z ned ) J 9) 2, C3 
if a 3.4 Go| & als - z 
40) Vip 8) 2. 25 
41 4 3 . |A af al A = | 
, ° M1 so). K. i A I Rs 
/ \ — Ae 2 Fe BeBe be * Be! 
, iit haes a 23 | Os 25 2. 
v¥, wll cd As Bebe Ag’ a By be al 
° A3 a 2 
| hey "lll g : ok ia? 
>) A 2. 
, <I A M0), , +[t2)e Aa, 
Fic. 9. Numbering of structural elements. 4 | pet | Ae 42 
16)| (lo, ,-/t2)- re, Se | 
Fol Fol | |A4 8 24 [ | 
i bid , , - a7| 73/+/t6) 
With the aid of these unit loads the self-equilibrating oa 
° m1: (3) Og_/ 2. 
systems L, can be represented in a tabular form. This ~*|“ os i 2 a, a 
will be shown in Tables 1-4. a Ned) Mit sa | os | Ome 
20) %, ffs) 4a v4 424 —_” z. 
° ° ° e - LL 
Fig. 9 represents a portion of a wing. Each panel, o | Zajfie) 2023 B23 | 
spar segment, and rib segment is marked with two Jf} ; rr. Oe me 
: +} (22)| (20), 1 -/20)p, | 2a & 
numbers, the first of which (f) corresponds to spanwise ‘ } ee Be | 
: z +] /23)| (2), ~/(2ty, 4Qp3 ¢ 
and the second (g) to chordwise numbering. Even CC } | “Ze 4 | 
e ° JI (24)| (7/ -(22)4/eZ | 
numbers are used for panels and odd numbers for spars _~ iad | ; wh 
P ° . e (Ps) — e a | 
and ribs. In the tables, the loading of a pair of panels, 5) My 4~ pes | 9% | the 
a rib, a pair of spar flanges, and a spar web is given in | 28) | Mey * oy | %*% | 
the form of factors by which the unit loads, noted at a a | 4-43 
the top of the column, must be multiplied. The cal- (28)) (ep y+ (elegy Ayohz | | 
culation of these multiplying factors can also be carried (ey) (25)ffe6) | a | / 
out in a tabular form. The numbers shown in paren- G0| (erffea) | fi | | 
theses after the factors refer to the numbers in Tables 5 3)| 9 -(0) | op-f> | 
and 6 where these factors can be found. Finally, the (32)| -(29)+(30) | s-ap/p | 
last column contains the numbers of the equations from (33) Gi}f+-/32) |  Op-Ap | 
° * | 4(7-G27e) | a 
which the factors can be calculated. 4)| (294/33) 
= | | 
Table 5 requires a little more explanation. The 44%7| 444 
factor R, for L; = 2.1 in Table 2 can be transformed as 4 4} —43)-4) 
follows: (a7)| 2(if(e6) | a aS | 
J 4) B 3-37) 
“) ay, 1 -( 2ashs3 ) ay ian ll — a | 1 ——- _ 
Az a3 Ayhg + azhy a3 
2ayhz ayhg sa ashy 
pee agg) ee is 
ayh3 + aghy ayhs + ashy 
Table 4 _ Web ir spor segment 2.3 
a a a aaa Se Table 6 ‘ er segments —_-/ 
Ly Nl Smo sp, rees| Mb Imb— ie cde lenceshar eh ég Zale eeeee oe ee 
et) -f he as | - (23) | t | ntereection pad / ? 5 | 
9 fob M13 + "33 - } : ae a I . 
12) ~ 3 | -Hoke |+ E8*"c23) + /K2),, _ wei A. | Aa | As Si 
fz a 7 | - ——-+— -—- — — — —4A - an see — 
32 Za | +()ee ~ “eg * "e203 ~ (o3)p3 - | (00) | (39), 1-139, | Kcmting: tegr 
2 c33 | wr « t t 1 ° 
23|+/ a, ol | (23) bot) | 3%. #89, 4 hes hes | or tl 
/3 * 3 2. t 1 in . 
Ae - = EB EEE | — fa2hpg | BE | +/O be | iti hal!) Mol [39 Mel hop | a clus 
24 i |. Mera *Pe33| . 43 2% | _ ig, oe 49! ll) fa%,, ie rae. | | f bes 
Ae33 7/23 Ze |“? nase 4 = [e} it o~—_d _ R 
2.5 = - | =f | (23) Gey) (al) | bores +hes) | a» 
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Table 7 Panel (fg) Table 8 Fb segment (fg) 
P OQ. Jo) ar | Ao t Fhe “be" ob “os *%4° Aes:| Fgs'| %s° “26° 
mm = —/- mm oe | 2 in  d ey -(22)5 Ma§ 1g" -le0)g|- 1g MS)8 
-(9) -(10), 1s @ | 
I og . 3) +e" -5' 02) -& | 7 |-67054 -4/i/ |3/ 024 07 
7 | +6 J 4" 4 
! “s | 
2 3 ] aa 
es +4 |43|\-2| —-f Load = 
Ca a | = we 
| +f ~— | Frac Fac’ Fea Fad (22) ond RIYJ resp vonish, if 
+3 ‘L235 (22); the spars Aaove ro break in 
Ee 
EE intersection nint Kg-l and 
Ye | +4 C ie on os? 
| M5 } fi. = Koel reso. 
42 i 
ae ST ee 
: : . r y \Ace | Fee 
and was quoted in this form in the table (see No. 7). Pew! 
The computation of this factor is then simplified, since = . 
it is observed that | fea), | +6" | 
. bat _| 
Ay A3 


(2 ]--(2-) = 
Az ay, Ao a3 ab 


(6) THE INTEGRATION OF THE PRODUCTS OF THE INTERNAL LOADS 





Two different load systems L,; and L; may cause forces in the same structural element, for instance, a rib seg- 
ment. The loading R,; of such a rib segment in load system L; can be expressed as 


R; — AiR i B,iRy + Ci.R, 


where A,,, etc., have the values listed in Table 2 in the columns R,, etc., in the row corresponding to L;. Similarly, 
the loading R,; in system L, is 


R; = A,;R, + B,,Ry + Cas 


Integration of the products of these two loadings yields 


| R.RdS = A,A,, / R,2dS + A,,B,; | R,RdS + AnCy | RaRAS + ByAr | RRS + 


B.By f R,*dS + Bucy f R,RadS + Cun f BRAS + CuBr f ReRedS + CyiCyy J R2dS (28) 


There are six different integrals in this expression which can be denoted by the following symbols: 


. 


/ RWS = KR. [ RRs = J Rras = Ry 
/ R,2dS = Ry, f R,RAS = rs RRdAS = Ree (29) 
/ RdS = Re. J R.R,dS f R,RAS = Rec 


Hence Eq. (28) can be written in the form: 
[Rras - A,;A,j;Raa + (A,,B,, is B,;A,;)Ra 2 B,:BryRo + (By iC; + C,By;) Roe a Cri CryRee + 
(Ar, vj + Cy1Ayy) Rae 


II 
II 
I} 


Since it is necessary to calculate all the definite in- and spar flange segments and one for panels and spar- 
tegrals °L,L,dS for all values of i and j, all products web segments. Some details of their calculation will 
of the coefficients listed in any one of Tables 1 to 4 in- be given in Section 8, while in the present section their 
clusive must be evaluated. These products must then values are assumed to be known. 
be multiplied by the corresponding unit integrals Ra, A convenient procedure for the systematic evaluation 
Ry, ete. Altogether, there are six unit integrals for rib of the products A,;A,;Raa, etc., is: 
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Table 9 Ent 2 Sper Tar _ 9) Table l2 Rearrangement of the integrals of products 
b= lbc* lab oe" “ag ae” | 36° Common spar segment raed (Fu even, g2007 ) 
) bse “foal, 4y 
L ‘BEY, ag Bey i aay page 36yf) ee te 
a; s Re oe oe OME es. 
1 f vA 7 “E> ge ES pape Oo | +f 7) 7) re) "O 
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| 7 ] 
(36); , . *2/ y | ? 
| 4 | 
fexy =9f = Ele | 
He = 0 -?| 47(1| S7(f S6(é 23(1| 224) (4(/| LIL | 
(38 )y f + 2 ue | | | 53(4| | | 
Po ah 2 | | £6(2| | 
(8g |-4 _,| 57(2| 5711) 46(6| s3(t| 13(2| 26(2| 24(4 
y, +f 67(2| 6711) 7él2| 43t} 4H16 F422) F6(/ 
BE)yo1 <Sf |+2 | | 7372 | | | | 
Gj / 2a ae | 2614 lost 25(4 oe 
(35lget = 5 o| ru 77u| Sela] 230| SHe| o4it $54) 22(f| 226) Wt) Ws 
(38yas oF [re 4 | } 4 5314] Tae | 52(4| + 4 
, TS |-7 | | | | 2 
rot x. l-4 47 } S4(2| S41) 4516) f2(f\ 12(2) 24(2 2rir| 
oe eos |+f 64(2| Salt) 7512) 4A) 426) 312) Flt} 
+3 | 2 
fe | -6' — art t ¢ | 13521 + + j 
02 | | | Fall) Fell) 65(4| FEL) FAS) 44(1) S41) 
+2 | | | 62(4| | 
} + t 4 } } } } -— + { 
a | | | | SH(2) 541} 
Table 10 Web in sbarsegment . (fg) £4 | | | | | | | | | 6£(2| 64(2} 
ra f ‘4 | Aaht 7] * . ae a a = ee a a a © LU) 7s 
MeS)y U4 Holy 4 | Maly foily [olye 1 
‘gpa 2a ge Pr grrr er are : ; 3 : oad 
me. My a Mah (49)q boély Wp. s- bey first column and in the first row, in the latter multiplied 
PEP Pg insite Sie rage ager ‘” : 
Mg bsiglole hale SF be Yor hDq by suitable unit integrals. Phe adjacent row and col- 
umn contain the numbering from 1 to 7. In the load 
"s~ far? | a <4 (ks Mob Mae | Ms’ systems containing several unit loads, case 6 is marked! 
Og MYg [hs yy ta Ag Mele Ase ys ‘ é gs : ; Ss, ~ shoes 
50)9} and cases a and ¢ are marked Case 4 of the ribs has 
i ° y ame | the ’ in the a case and the "in the c case. It might be 
— A i I Nie LR te Sace SonA M noted that the only case containing three unit loads 
la ° . ° 
(sig f namely, case 4 of the ribs—contains all three of 
Mob «2 , ( 
f lg Zs them only when both spars have a break at the rib 
47, *3 = 
: 4 concerned. 
b5)g oad The quantities of the first row and of the first column 
40] are now multiplied, anc e products are denoted by 
+/50) Itiplied, and the products are denoted by 
* ees oe 
bA)y | ind double symbols containing the two original symbols. 
4 , o . . ‘ a] 
a “6 The order of the symbols is unimportant, so that 46 
; i 
0 - ; : cei 
209 | . means the same as 64. As far as the rib is concerned 
one can write, for instance, . 
Taxblo Mf fearvengemant of the integrate of premicte (46) = (64) = 6!4! — 6’4! + 614’ — 6/4’ — 614” + 64" 
Commer 4 segment Fy ~ bts ~ 00d, Ju eve 
Lg =e a we Sa > STAR] The algebraic signs must be preserved in the equations. 
al w a Ww -iw | r T ar - "f : 
ah eer rr ar oe er oT 1 SHE SCOT he above sum of products represents JS R:R,dS for rib 
| 497, ¢ 4 jes | Lvs | +. / oO of / 5 ee ° Q« 
4-3 | ele) 76 7 as iar cere segment 3.2 and for the load systems L; = 3.2 and L 
2 | pete | B67 | | — 1.3. 
0-2 3E(4 | F6(/ 25(%| 25(1 ° . 
5 janet — I | 2816 | 334 It should be noted that the calculations outlined 
(214 b 2 46 . * 
2 alae aed eal (9 OE ONE asiolasu lor |re(2\ seu 2r(2} above must be carried out separately for each structural 
*2 + ; + ; + St = : y 
: Fah, ial iaiad element such as rib segment, spar segment, etc. Each 
2 47 (- 6 14( 23I(4 Ot fo(32 P ° - . ° - 
eo “le COU 44 (1 | 33(4) 3310 | guid | S5(4) S514) 11(4 | 2244/22/41 table contains, therefore, the integrals of the products 
Stee ‘“) for a given structural element for different combinations 
f 7|67(2 | 7602) 76(1 | 3401 | 43(2 | 4314 | 2416 ° 
hl aes | | Bale | 7512 | 754! a4U4 | 212 | walt | 3 of the load systems L; and L;. On the other hand, the 
“3 | | 6/(2 . '< ~ . ‘ ° 
Am? 53(4| 531 coefficients f/ L,L,dS in Eq. (3) necessitate the deter- 
02 3(4| 63(f $2(%| $2lt ; é i. : . bis: 
+2 | 62(4| 62/1 | mination of the integrals of the products for dz/ferent 
“t43 73H(2) 7h x ss } : ‘ ‘ ; : : 
aaa 1 ___izel2|7e/| | structural elements for a given combination of the load 








The load systems are numbered from 1 to 7 in the 
order in which they are listed in Tables 1 to 4 inclusive 
(in the case of the panels, they are numbered from 1 to 4 
inclusive). Since most of the coefficients belong to the 
unit loads 6 and since most of the coefficients belonging 
to the unit loads a and ¢ are equal to unity, the compu- 
tation of the triple products is not too laborious. This 
Tables 7 to 10 inclusive. The coef- 
lables 5 and 6 are listed in the 


‘an be seen in 
ficients determined in “ 


the values contained in 
as to 


and L;. 
lables 7 to 10 inclusive 
permit the summation of the values corresponding to 4 
and L,; over 
Tables 11 


systems L Hence, 


must be rearranged so 
given combination of the load systems L, 
all the structural elements. This is done in 
and 12. 

An example will clarify the procedure. The definite 
integral {L,L,dS is evaluated for L; = Ls.4 and L, = 
Lo.5 5(odd) and g = 4 (even). 
Hence, Moreover, Af 


Consequently, f = 


Table 11 must be used. 
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»9—5 = —3and Ag = 5 —4= +1. Inthe relevant 
row only one number is listed: 37(2) under R.—2.0. 
Consequently, the value under 37 in the table for rib 
Rs-2440 = R34 must be filled in. The number in 
parentheses indicates that two values must be read off 
the tables and added. They are obviously the values 
under —3!7 and —3’7, both with inverse sign. When 
the row of L,; contains several numbers, then the corre- 
sponding values must all be summed up. It is clear 
that interchanging L; and L, does not change the result. 
In the example, Af becomes +3 and Ag becomes —1. 
In Table 12 one finds again 73(2) under R. +1. — 1. 
The corresponding rib is Ro;5-1 = R34, as before. 

It can be seen, therefore, that it is sufficient to calcu- 
late for each equation the coefficients situated at one 
side of the principal diagonal, including those on the 
principal diagonal, in the matrix form given in Table 
I. If each of the coefficients is determined, then 
those to the right of the principal diagonal can be used 
as a check of those to the left of it. In the case of the 
coefficients on the principal diagonal, the number of 
values to be added is always 1, 4, or 9. For instance, 
in the case of the spar flanges the combination 66 means: 
6!6! — 616’ — 6’6! + 6’6’, which is the same as 6!6! — 
2(616’) + 6'6’. 

Moreover, Tables 11 and 12 contain in parentheses 
the number of values that need to be taken into account 
for the ribs in the general case of broken spars. If these 
hints are observed, no numbers will be overlooked. 


(7) Tue EXTERNAL Loaps 


The internal load system L» is an arbitrary system that 
balances the loads applied to the structure. If this is 
kept in mind, there are basically two possible choices 
for a suitable load system: 

(a) From an estimate of the statically indeterminate 
quantities, the load system JL» is established in such a 
manner as to obtain the best possible approximation to 
the actual system L. The advantage of this procedure 
is that the structural elements can be proportioned in a 
satisfactory manner from the loads Lp alone, without 
making use of the loads L;. Moreover, the statically 
indeterminate quantities are all relatively small, and 
for this reason no great accuracy is required in their 
calculation. 

(b) The load system Lp is established in such a man- 
ner as to involve the least possible number of structural 
elements. The advantage of this procedure is a ma- 
terial saving in the work of computation when the defi- 
nite integrals S Loh dS are evaluated. 

Which of the two procedures is preferable has not 
yet been established. Moreover, the choice may de- 
pend on the structural arrangement, as well as on the 
individual experience of the stress analyst. 

The calculation of the loads Lo is a statically deter- 
minate problem. It can be solved without difficulty 
with the aid of the formulas given in Section 4 and the 


elementary laws of statics. The evaluation of the in- 
tegrals fLoL,dS differs from that of the definite in- 
tegrals discussed earlier insofar as the loading Lp cannot 
always be decomposed into a limited number of unit 
loads. In particular, loads applied between the inter- 
section points of spars and ribs cannot be reduced to 
one of the unit loads of Section 5. This difficulty can 
be overcome in the following manner: 

The loading Ly is decomposed into two parts: (1) 
The loading Zo’ that would prevail if all the intersection 
lines of ribs and spars were rigidly fixed and the external 
loads were applied to this system; (2) the loading Lo* 
arising from external loads, equal and opposite to the 
reaction forces and moments of loading Lo’. This load 
system is balanced by the actual reaction forces and 
moments. 

Shear flows in the panels are obviously nonexistent 
in load system Lo’, since the corners of the panels are 
not displaced. Each segment of spar and rib is a beam 
with two rigidly fixed ends whose internal loads depend 
only upon the external loads applied to it. 

In load system L* the external loads are all applied 
at the intersection points of spars and ribs. This system 
can therefore be decomposed into the unit loads of Sec- 
tion 5. Consequently, the definite integrals f-Lo*LjdS 
can be evaluated according to the procedure discussed 
earlier. 

If the statically indeterminate quantities correspond- 
ing to Lo* are calculated, the load system L* can be ex- 
pressed in agreement with Eq. (1) as follows: 


L* =Lo* + O X-AL, 
¢=] 


On the other hand, the load systems L’ and Lo’ are 
identical because each spar and rib segment can be con- 
sidered as an independent structural element. Hence, 


L=L*+L’ =1Lh* +l) + > X*L, = 
i=1 
lo+ > XL, (30) 
i=] 


The computation of the internal loads corresponding 
to Ly’ necessitates a statically indeterminate calculation 
of each spar and rib segment. On the other hand, they 
amount only to a small percentage of those correspond- 
ing to Lo because the same loads act upon much smaller 
units in the former system than in the latter. For 
this reason an estimate of the fixed end moments with 
an accuracy of 10 to 15 per cent causes an error amount- 
ing to an extremely small percentage when referred to 
load system Lo. Moreover, an error in this estimate 
has no effect upon the loading Lo = Lo’ + Lo*, since it 
cancels out when the two component systems are added 
together. An influence is present only in the definite 
integrals /"Lo*L,dS. Since great accuracy is not re- 
quired in the evaluation of these integrals, as already 
mentioned in Section 2, an estimate rather than an 
actual calculation of the fixed end moments should be 
satisfactory for practical purposes. 
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External loads on part of rib segment. 


Fic. 11. 


(S) DETERMINATION OF THE UNIT LOADS AND THEIR 
INTEGRATION 


In Section 6 the integrals of the unit loads were as- 
sumed to be known. Now it will be indicated how these 
integrals can be calculated. At the same time some re- 
marks will be made regarding the internal loads in the 
ribs. 

Structural details of the various elements vary con- 
Hence, it is impossible to give general rules 
The suggestions that follow permit 


siderably. 
for the integration. 
an integration without undue difficulty. 

Considering the statements in Section 2 concerning 
the required accuracy in evaluating the load integrals, 
the following assumptions are made for the purposes of 
integration: (a) The thickness of the sheet covering is 
constant for each panel. (b) Variations in flanges and 
webs are assumed to be continuous. 
in the cross-sectional properties are replaced by esti- 
mated equivalent continuous variations. 

The unit loads shown in Figs. 5 to 8 inclusive will now 
be considered individually : 


Discontinuities 


(a) The Panels 


There is only one relevant unit load case P,, in ac- 
cordance with Eq. (11): 


dp = Gm2((1 — a*) (1 -_ doa) *] 
where, according to Eqs. (24) and (9), 
dm = 1, (2A m2) = | (ashy + ash) (31) 


The value of the sum a,/3 + a3/; is to be taken from 
row 6 of Table 5. 
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With the aid of Eqs. (29), (2), and (11) one obtains, 


for a pair of panels, 


?, 2206 oro 
/ P,2dS = f Se a = 
: Gte 


9 Ym2"(1 = Vets fi dry ie 
Gte -1 Q— A2ae)3 . 


The integration of this expression, with due consid 
eration of Eq. (31), yields 


: = 2asby [G(aihs = 3h) *te} (32 


A somewhat better approximation is reached by tak 
ing for dab. in Eq. (32) the area of the panels between 
the rivets connecting them to the adjacent spars and 


ribs. 


) om = 2A 42/[G(arhs — a3h;)*te] (de 


(b) The ribs 

The internal loads in the ribs, caused by the external 
loads shown in Fig. 6, depend greatly upon the shape 
of the rib. This is particularly true with regard to the 
shear force in the web which greatly depends upon the 
difference in slope of the upper and lower flanges of the 


rib. 


The slopes cannot always be determined accurately 
enough from the drawings, particularly not when they 
are not drawn exactly to scale. It appears advisable, 
therefore, to fix the shape of the rib as a polynomial of 
the third order of the chordwise coordinate. This 
should be done both for the total height between the 
cover sheets of the wing and for the effective height /, 
measured between the centroids of the flanges. Each 
function can be specified by means of four ordinates, as 


indicated in Fig. 10. 
The bending moments and the slopes can then be 
given in the following form: 
~ 1 > 
O(hy + 3he oh 3h3 + hg) m=1 


4 


M, Ri my m (34 


w 
or 


hep’ = (1/22) D> Rempltem (: 
m=1 


where m is the station number for the ordinates and » 
is the station number for the cross section under con- 
sideration. The coefficients Rimp and kemp are fixed 
integers, presented in Table 13. 





Table 43 

Factor kim k2m.p 

pol 1|2)}2|4] 7] 2)/3 14 
#9 1427\427)\4 GIi-Lf |\4#18|-9)|+2 
0 l+ @le32i¢ B|—- 2|/- 3/+ 6|-7 
— 5|+79|+ 9|+7|-6|+3|42| 
oO ra] O O|- 2|+ 9|-1al +i] 
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With the aid of Fig. 11 one obtains, finally, in cases R, and R-.: 
Sp = gr(2hy — hey) + My (h' cp/Itep) 


Substitution of 


gr = 1/(2A,) = 4/[a(i + 3h2 + 3h2 + hy)] 


results in 


4 
aid a(hy + Sho + 3h3 + hy) 


In case R, one obtains 


4 
S, = (2hy — he) — 


a(hy + 3h2 + 3h; + hy) 


The integration can now be carried out. 


f= dx = [ 
El a 


Ga 


f GaGa dx = f 
Gt 


and 


= (2h, — 


hy + 3h2 + Shs + Ig 


4 < 
(= binabn)( hs cvten) (36) 


The 


hep) + 


4 4 
( } Rip vin) :. benaltem (37) 
+ = mot ' 


4 Thep 


One has to calculate 


MaMy _ 


~ eo -dx = Rj, ete. 
2) EA sh? 


SaSplt 
= _—< dx = Ra»; etc. 
Gth? 


where A, is the cross-sectional area of the flange, ¢ is the thickness of the web, and hg is the height of the web be- 


tween the connecting rivets. 


On the assumption that the integrands can also be expressed approximately by 


polynomials of the third order, one obtains for Ry» and Ra»: 





iT ALM, M,M M,M, M,M. 
Rai = - {| | + 3 | ee | +3 ia: a 4 + | —* | t (38) 
8 1L(1/2)EAsh2 J; (1/2) EA jh? 2 (1/2)EA ,jh,? Js (1/2)EA,h,?\ 
‘SaSoh S,Syh S, Soh S,Syh 
R, —- a | (> b ‘) ' 3 ( ab ‘) 3 (3 b ‘) t ( ab ‘| (39) 
== 3b Game), + FN Gant). t 2 Gut ls > ome) 


where \/,, M,, S,, and S, are known from Eqs. (34), (36), and (37). 


c) The Spar Flanges 


Internal loads caused in the spars by bending are 
given completely in Fig. 7. It should not be forgotten 
that the definite integrals involved represent the largest 
contributions to the coefficients in Eqs. (6). The ac- 
curacy required for the unit integrals F is therefore 
greater than that necessary for the other definite inte- 
grals. 

Asimplification results from the absence of any param- 
eters in the load systems F,,, F;, and F, (see Fig. 7). 
The definite integrals can, therefore, be given immedi- 
ately in the form: 


. MM, l (; ) = Ravm 
Fo = ~ i= .- = (40) 
’ J ae ”6| CUR GS S: T 


Here ” is the number of equal subdivisions of the 
Spar segment. For a given value of n, the multiplying 
factors kgm can be calculated once for all, just like the 
factors Rimp and kom, of the rib segments. They are not 


presented in order to save space. This procedure is 
again based on the assumption that the variation of the 
bending rigidity can be represented by a polynomial of 
the nth degree. 

Any accuracy desired can be obtained through a suit- 
able choice of m. Probably x = 2 will mostly be accu- 
rateenough. Insome calculations, however, » = 1 may 
suffice. When the bending rigidity is constant, one 
obtains with m = 0 the well-known expressions: 


il (i) 
3EI \b 
ra 2 (t), ra 8 (1) 

6EI \b 15EI \b 


The multiplying factor (//b)? is due to the fact that 
in the unit case the X components of the bending 
moment are put equal to unity, and thus the actual 
moments have the value of //d. 


Fes = Fay = Fre = F ee = 


(41) 
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(d) The Spar Webs 


In cases W, and IW’, the loading follows from Eq, (15): 


= l l —~ 2 = 
— 2/ —_ ha)? 
al l — a’ _ aden =s\t (49 
ay —6L(1— Xa)? (y¥ — rd)2]) . 
where @ and 8 represent the values corresponding to 
the cell immediately ahead in case IV, and the cell im 
mediately following in case IV. 
In case W, the load can be calculated from Eq. (11 
after substitution of 6, y for a: 


Quo = (1/2An)[(y? — &)/(y — XO)? ] (43 


Eq. (42) can be rewritten in the form: 


lj v— 6 
‘ 2A \(y — Nb 
a(l — y) — (8 6 | — a y? — 6 }) 
_ = = olf | t4 
oT = ¢ (1 — Aq) (y — XO 


where the first term in the braces is identical with the 
corresponding expression in Eq. (43). Substitution of 
\ = —1, 0, and +1 in Eqs. (43) and (44) vields, after 


some manipulations, 


Quoi = (1 2A »)[(4 
2A ») [1 — (6/y¥)° 15 


Qun2 o (1 
(+y — 6 | 


Qua = (1 9A.) I (4 1 §) y 


6 es is 6)]] 


a 
x 
~ 
Po | 
| 
v7 
—— 
> « 
, ~ 
o 
—— 
— 


~ 


2A 


l vy + 6 ws 2[a(l — vy) — (8 — 9)}I 
» io (1 — a)(y — 6 y 


It is easily seen that the second terms in the brackets 
of Eqs. (46) are extremely small as compared with the 
first terms. Furthermore, they change sign for —1 — 


A\— + 1. In evaluating the integrals 


Wa = / (qa9rita/Gt)dl, ete. 


they might be neglected. Only the integral II’,, there - 
fore needs to be determined. 


the panel integral P,,,,,, 


It yields, in analogy to 


Woo _ A wa (2A ,)?Gt (47) 


where Aya is the real area of the spar web between the 
rivets connecting it to the spar flanges. 
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1949 


(9) THe ULTIMATE LOAD 


The statements in Section 2 concerning the accuracy 
of the statically indeterminate calculation yield that a 
computation, based on the Theory of Elasticity, does 
not give satisfactory results in the case of the ultimat 
load, since at that load the modulus of elasticity changes 
rapidly and differs materially in the various structural 
elements. A more accurate stress distribution can be 
obtained, therefore, if on the basis of the first calcula 
tion, the rigidities of the structural elements subjected 
to the highest stresses are reduced. Such a procedure 
is not convenient, since it requires a new solution of the 
system of equations for each case of loading and after 
each approximation of the real stress distribution 
Each new system of equations has a new set of coef 
ficients. 

A better solution is available in the method proposed 
by G. C. Best in reference 2. The principle underlying 
Best's solution can be used in the present problem as 
follows: 


Let us assume that at a certain station the stresses 
calculated cause a reduction in the value of the modulus 
This reduction causes an 
The same in 


of elasticity from £ to nF. 
increase in the strain by a factor 1/y. 
crease in the strain can be achieved through an increas¢ 
by a factor 1/n of the internal loads, while the modulus 
of elasticity and consequently the rigidity are kept un 
changed. , 

As a result, the internal load L must be increased 
locally by the amount L[(1/n) — 1] = L[(1 — 9)/n] in 
order to obtain L(1/n). Of course, the value of 7 varies 
with structural element and loading. 


This increase can be taken care of by increasing Ly t 
Lo + L[(1 — )/n]. Hence, the statically indetermi 
nate calculation must be repeated with the same coef 
ficients $L,L,dS but with the new set of right-hand 
members /°[L(1 — n)/n]L)dS instead of f-LoLdS 
The redundant quantities obtained in such a manner 
are then used to correct the loads L determined 


earlier. 


In general, these calculations will result in a reduction 
in the stresses of the most highly stressed elements 
It is advisable, therefore, to estimate 7 too high in the 
calculations. Then the values are likely to be more 
nearly correct for the stresses obtained from the second 
calculation. The procedure described is basically an 
iteration procedure, which, however, is not too laborious 
because the matrix of the equations remains unchanged 
and only the right-hand members are replaced by new 
ones in each calculation. 


It must not be forgotten that a structural element, 
lightly loaded after the first calculation, may become 
heavily loaded after the corrections. For this reason 
the correction has to be carried out for all the elements 
of the structure. 
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APPENDIX 


Tables I and II show the solution of a system of four 
simultaneous linear equations according to Gauss’ 
method. Table I does not require much explanation. 
It results in Eq. (4), allowing the solution of Y,. With 
3) the unknown YX; can then be determined; with 
Gauss’ method was thus presented. ! 


Eq. (3 
Eq. (2), Xo; ete. 
However, it is usually convenient to carry out the 
calculations with the right-hand members in a separate 
table, especially if several sets of right-hand members 
are to be dealt with. This is shown in Table II, in 
which the multiplying factors of Table I are written in 
matrix form with the right-hand members below. The 
first element is repeated and (following the arrow) mul- 
tiplied successively by the matrix coefficients of the 
first row. Then the second column is added, and the 
sum is multiplied successively by the matrix coefficients 
of the second row, etc. These computations are identi- 
cal with those of the last column of Table I. 

The determination of the unknowns leads to the fol- 
The first coefficients of Eqs. (1), 
and the 


lowing computations: 
(2), ete., are written down in the next row, 
sums calculated are divided by the coefficients below. 
Next, the last element is repeated and multiplied suc- 
cessively by the matrix coefficients of the row below. 
Then the next to the last column is added, and the sum 
is multiplied by the coefficients of the next to the last 
row, etc. The sums so obtained are the unknowns re- 
quired. 

An investigation of the time required for these calcu- 
lations gave the following results: 

If one assumes that 54 sec. are needed for positioning 


the slide rule, 30 sec. for reading it, and 6m sec. for 
adding up (m + 1) members, including the time for 
writing down the numbers, then the following figures 
are obtained: 


Time required for Table I: 


T, = n(n? + 2n + 3)/5 min. 


Time required for Table II: 
T, = (3n? + 13n + 9)/5 min. 
Time required for a general check: 
T3 = 


3n(n + 4)/5 min. 


In the time requirement for Table I, the time neces- 
sary for the last column is not included because the cor- 
responding calculations are carried out in Table II. 

The basic time requirements for the fundamental 
manipulations (54, 30, and 6 sec.) may appear too high, 
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£guations: Table I 
Pi :+S$987%,-(523kpe IPR; =— 22 
8:- ~ 152.3%, #3775 Xp = “5.125 + (20%, =+608 
Cre 92% = 45.1 Xp + 457.02 = F049, = —428 
D: + 120% -IOMGR + IWIX, =—- 224 
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but they seem to be necessary in reality, The resulting 
time requirements were found correct, so far as order of 
magnitude was concerned, in a number of practical 
calculations. It should also be noted that the compu- 
tations can usually be carried out accurately enough 
with the slide rule, even when -the number of the un- 
knowns is large, because many coefficients become ex- 
tremely small or even vanish. 








On a Complete Solution of the One-Dimen.- 
sional Flow Equations of a Viscous, Heat- 
Conducting, Compressible Gas’ 


MORRIS MORDUCHOW? anv PAUL A. LIBBYt 
Polytechnic Institute of Brooklyn 


SUMMARY 


The hydrodynamic equations describing the steady, one- 
dimensional flow of a viscous, heat-conducting, compressible 
gas are treated. The equation of state for a perfect gas is as- 
sumed valid, and the coefficients of specific heat are assumed 
constant. A complete integral of the energy equation is found 
for a Prandtl Number of */,, which is a good approximation for 
the actual value of many gases over wide temperature ranges. 
A particular integral contained in this complete integral is (u?/2) 
+ ¢cpT = aconstant, which is a frequently used one-dimensional 
energy equation. The restrictions on its use are rigorously de- 
rived and stated. With the complete integral of the energy 
equation, three different types of solutions which depend on the 
boundary conditions are shown to be obtainable. The usual 
shock wave solution (satisfying the Rankine-Hugoniot relations) 
is only one of these three types. 

Several features of the shock wave solution are shown. For 
example, the entropy has a maximum value at the inflection point 
in the velocity distribution. The validity of the continuum 
hypothesis is examined, and it is shown that the thickness of the 
wave as defined by Prandtl becomes infinite as the initial Mach 
Number of the flow becomes infinite provided that'n, the expo- 
nent in the temperature-viscosity relation u/uo = (T/T)”, is 
greater than '/>, Furthermore, it is concluded that the contin- 
uum theory may give reasonably correct results for flows with 
Mach Numbers up to 1.38. 

The other two cases of flows contained in the equations are 
terminating flows—i.e., flows that can be valid only in a finite 
region. Their mathematical nature is demonstrated by two 
numerical examples, which indicate compression shocks followed 
by expansion waves. However, the physical significance, if any 
exists, of such solutions remains open to question. 


INTRODUCTION 


Siw PURPOSE OF THE PRESENT PAPER is to give a 
concise, and yet general, theory of the 6ne-dimen- 
sional flow of a real continuous fluid in order to derive 
all of the mathematical implications inherent in the 
equations governing such flows. Solutions with uni- 
form and nonuniform conditions behind a shock front 
and the relationship between these two types of solu- 
tions are explicitly shown. 

As an approach to the mathematical description of 
the structure of a shock wave, an analysis may first 
be made of the one-dimensional steady flow of a con- 


tinuous medium. Such analyses have already been 


Received February 10, 1949. 
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conducted by several investigators, but these have, in 
one way or another, been limited in generality. In 
early investigations,’ for example, the flow of a (hypo- 
thetical) viscous fluid without heat conduction, and of a 
(hypothetical) nonviscous fluid with heat conduction 
were treated. Taylor and Maccoll* later gave an ap- 
proximate solution for the flow of a fluid with both vis- 
cosity and heat conduction, valid only for a weak shock 
wave. In 1922, Becker*® obtained an exact solution of 
the one-dimensional equations of a real fluid with a 
Prandtl Number of */;, by assuming the expression for 
the temperature to be a quadratic in the velocity and 
then determining the coefficients by substitution into 
the differential equations. Thomas‘ afterwards briefly 
extended Becker's investigation by using variable coef- 
ficients of viscosity (u) and of heat conduction (), as- 
suming uw and & to be proportional to 7"/,. The latter 
investigation, which introduced directly the mean free 
path, showed the increase of shockwave thickness due 
to variable values of » and k, although, as shown in the 
present paper, some of the mathematical conclusions 
that Thomas reached must be modified if a different 
law for the variation of u and k with 7 is assumed. 


All of the aforementioned investigations had been 
restricted to the case of uniform-flow conditions both 
ahead (x — —o) and behind (x — +o) a shock 
front. Reissner and Meyerhoff* have recently pointed 
out, however, by means of a small-perturbation pro- 
cedure, the existence of solutions to the shock wave 
equations which do not represent uniform conditions 
behind the shock front. 


In this paper the one-dimensional flow equations are 
solved rigorously and completely for a constant Prandtl 
Number of */;. The existence of solutions correspond- 
ing to flows with nonuniform conditions behind the 
wave is clearly proved and is demonstrated by two nu- 
merical examples. The thickness of a shock wave in 
terms of the mean free path of the gas molecules 1s 


considered. 


SYMBOLS 


0 = subscript denoting values at x —~ — © 

1 = subscript denoting values at either x = x; (or x = 0) 
oratx—> + o 

C1, C2, Cs, Ge. Ge. C,” = 

Cp, Ce = specific heats at constant pressure and at constant 
volume, respectively 


constants of integration 
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FLOW EQUATIONS 

H = the quantity (u?/2) + c,T 

k = coefficient of heat conductivity 

M = Mach Number 

n = exponent in the viscosity-temperature relation 4 = 
uo( T/T)” 

p = pressure 

R = gasconstant R = cp — & 

§ = entropy 

i = absolute temperature 

ty = thickness of a shock wave, according to Prandtl defi- 
nition 

ty = value of ¢; for n = 0 

ti = thickness of a shock wave, according to Taylor- 
Maccoll definition for » = 0 

i" = thickness of a shock wave (Prandtl definition) for 
n ~ 0 

tr = thickness of a shock wave in units of local mean free 
paths. 

u = velocity of flow (in x direction) 

V = dimensionless velocity u/1o 

x = distance coordinate in direction of flow 

: = exp [2 cyMit/(y + 1)] 

a = parameter defined by (y — 1)/(y + 1) + 2/(y+ 1)X 
M,? 

6 = aconstant, 0.35 

1 = ratio of specific heats cp/Ccp 

K = constant defined by Eq. (11b) and having the value 
1.36 for air (y = 1.4) 

d = mean free path of the gas molecules 

mi = coefficient of viscosity 

p = mass density of the fluid 

1 = dimensionless temperature 7/7 (at x = x) 

: = dimensionless distance coordinate x/Xo 


EQUATIONS FOR ONE-DIMENSIONAL FLOW 


The steady, one-dimensional flow, parallel to the x- 
axis, of a viscous, heat-conducting, compressible fluid 
is described by the following hydrodynamic equations, 
which, in the kinetic theory of nonuniform gases, form 
a second approximation to the Boltzmann equation: 

Conservation of momentum: 


du dp 4d “*) 
: Easechagaed -)}= la 
ia dx + dx nl! dx ° (la) 
Conservation of mass: 
(d/dx)(pu) = 0 (1b) 


Conservation of energy: 


pu + pS — au") ~2(#Z) = 0 C00) 


dx ax 3 dx 7 dx dx 


Equation of state: 


p = ple, T) (14) 

where 

p = mass density 

u = velocity 

p = hydrostatic pressure 

& = coefficient of viscosity 

E = internal energy 

T = absolute temperature 

k = coefficient of heat conductivity 
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For the usual ranges of pressure, density, and tem- 
perature encountered in fluid flows, the coefficients of 
viscosity and heat conductivity may be assumed to be 
functions only of the temperature. The internal 
energy will be a function, in general, of p, p, 7. The 
coefficients of specific heat at constant volume and 
pressure (c, and c,, respectively) are considered here as 
constants.* The four equations [Eqs. (la)—(1d)] de- 
termine the flow variables u, p, p, and 7 as functions 
of x. 

For mathematical simplicity the particular equation 
of state which has been used here is that for a perfect 
gas—namely, 


p = pRT (le) 


where R = the gas constant, R = c, — c¢,. With this 


equation of state, 


dE = cdi (1f) 


The equations can be integrated as follows. The flow 
is assumed to be in the positive x direction, and three of 
the five constants of integration will be determined by 
conditions at minus infinity, which will be denoted by a 
subscript zero. Eq. (1b) can be integrated directly to 
give 

pu = polo (2) 
Eq. (2) may be substituted in Eq. (le) to give 
p as poloRT ‘u 


or upon differentiation 





a 
dx u dx 


(3b) 
udx 


With Eq. (3b), Eq. (la) gives 
p du P du . pet 4 u = *) (4) 


—-——(s- 
3 Pollo dx dx 


pollo dx dx dx 


Furthermore, Eqs. (2) and (4), with Eq. (lf) and the 
relation c, = c, + R, can then be used to transform 


Eq. (1c) into 
du dT 4 dfd *) (*) | 
ax t % dx =| aC dx Th dx 


l ( ~~) 
——(k—)=0 
pollo dx dx 


du dT 4 “A “*) - d (eT) a 
m dx il dx  3pouodx — dx polo dx\ dx 


(5b) 


(5a) 





By a single quadrature, Eq. (5b) gives: 
* Reference (6) indicates a change of 16 per cent in Cc, for air 
for a range of temperatures from 0° to 1,600°F. 
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u* k (4uc, du dT ; 
-+¢,7 — - a: “— c = €, (5c) 
2 . Cppotto\, 3k dx ? dx 


where C, is an arbitrary constant. 
If the Prandtl Number 


ucp/k = */4 (5d) 


then Eq. (5c) becomes 


k d/([(u? 
(“ +a) =C,  (5e) 


Cppotto dx \ 2 


P 
u* = 
= + oI — 
Eq. (5e) can be integrated to give 


u* Be a ; dx 
ry + fal = Cy + Co Exp | Cypollo ry (6a) 


Eq. (6a) may be regarded as a complete integral of the 
one-dimensional energy equation for a fluid with a 
Prandtl Number equal to */,.* A necessary (though, 
as will be seen later, not sufficient) condition that the 
flow be uniform at plus infinity is that C, be equal to 
In that case Eq. (6a) reduces to 


(u?/2) +oT = CG 


zero. 
(6b) 


Thus it is seen that the well-known Eq. (6b) is an 
exact and complete integral of the general one-dimen- 
sional energy equation of a viscous, heat-conducting 
fluid under the following two conditions: (a) the 
Prandtl Number is */;; | (b) the flow is uniform at plus 
as well as at minus infinity. Condition (a) is necessary 
so that Eq. (6b) be exact, while condition (b) is neces- 
sary so that Eq. (6b) be complete. Although Becker* 
essentially obtained Eq. (6b) for a Prandtl Number 
equal to */, by assuming 7’ as a general quadratic func- 
tion of uw and by determining the coefficients in this 
function so as to satisfy the governing differential 
equations, it appears not to have been generally recog- 
nized that this quadratic was the familiar Eq. (6b).7 
Moreover, by the assumption of the temperature as a 
quadratic function of velocity, Becker's solution was 
implicitly restricted throughout to the case of uniform 
conditions at plus infinity. Therefore, his paper does 
not contain Eq. (6a)—the complete integral of the 
general one-dimensional energy equation. 

It must be noted that condition (a), requiring P, = 
3/4, is based mathematically on the factor */; appearing 
in Eq. (5c) and, therefore, on the Navier-Stokes mo- 
mentum Eq. (la) and the energy Eq. (lc). Hence, 
condition (a) is based on the physical assumptions used 
the Navier-Stokes equations. 
(i) The components of the stress tensor 


in the derivation of 
These include: 
are related linearly to the components of the strain 
tensor: and (ii) the hydrostatic pressure is equal to the 
arithmetic average of the normal stresses. 

* It may be pointed out here that Eq. (6a) is valid for variable 
values of u and & provided only that the Prandtl Number remains 
constant and equal to */;, k being an implicit function of x 
[ke = R(T) = R[T (x)] = ki(x)]. 

t This may have been due to the somewhat obscure notation 
used by Becker. 


SCIENCES-NOVEMBER, 1949 

The case of uniform conditions at plus and minus 
infinity appears to be the only one considered prior to 
the work of Reissner and Meyerhoff,> who indicated 
the possibility of solutions with nonuniform conditions 
at plus infinity. With Eq. (6a) a first order, nonlinear, 
differential equation for u as a function of x can be ob- 
tained; solutions of this equation with C, negative, 
zero, and positive will be discussed. 

If Eq. (2) is substituted into Eq. (la), the latter may 
be integrated once to give 


Pollou + = (4/3) u(du dx) = C; (7a 
where C; is an arbitrary constant. Eggs. (3a) and (6a), 
combined with Eq. (7a), with y = c,/c,, lead to the 
following first-order differential equation for 1: 
4pup du <ytil, C3 
<n ee gt “u = 
3 pollo dx 2y pollo 


— |] 1x +1 
d ines C2 exp seus f | + z nes C; (7b) 
Y k Y 


C, and C; can be determined from given conditions at 
minus infinity, where all the flow variables are assumed 
to have finite values—that is, wu = uo, T = To, p = pr, 
p = po, and d/dx = 0. Eqs. (6a) and (7a) then yield 


Ci = (140? /D) aa Cyl o (Sa 


and 
C3 = Polo” + Po => po(Uo” oa RT») (8b 


In treating Eq. (7b) it was found convenient to non- 
dimensionalize u with respect to uo and x with respect to 
Xo, the mean free path of the gas molecules at minus 


infinity. Thus, let 


VY = — = x/Xo (9a 


u/Uo; 


Moreover, certain relations from the kinetic theory of 
gases may be used; thus, one can write® 


where 
6 = aconstant = 0.35] 
¢ = mean molecular velocity 
= VSRT/x 
\ = mean free path of the molecules 


Furthermore, both theory and experiment indicate 
that u = pw(7); in fact, it is usually satisfactory to as- 
sume for gases that 


= po( T° Ze)" Ye 


where n is a positive constant depending only on the 

t Theoretical values of 6 vary from '/; to 0.499, depending 
the nature of the analysis. The value used here is that of Tait 
The effect of using a different value for 6 would be to change the 
scale of — (and therefore the shock-wave thickness in terms of 
mean free paths) in direct proportion. 
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With Eqs. (8a), (Sb), and (9a)—(9c), 


gas in question. 
Eq. (7b) becomes 


*\n V 
(=) = — «Mo (V — 1)(V — a) = 


T; dé 
y¥ —1\C 2xy Mo Mo Ze\" 
9, 5 EE a 10 
2a > : ye exp | 4 r dé| (10) 
where 
Vf) = Mach Number at minus infinity 
= Uo WV yRT» 
ter = 7— + —_— a 
a = aparameter = —— —— a 
i parameter yl “et 7 Me ~ 
« = aconstant* = + 1) yz (11b) 


FLOW WITH UNIFORM CONDITIONS AT PLUS INFINITY 


General Solutions 


With C2 
flow be uniform at plus and minus infinity, Eq. (10) 
can be integrated directly to give 


I i (T/To)"VdV ; 
t= — : ——+Q(, 
«Mo (J = 1)(1 sary 


It may be pointed 


= 0, which is a necessary condition that the 


(12a) 


where C; is an arbitrary constant. 
out that one obvious solution satisfying Eq. (10) with 
C, = 0 and the boundary conditions is V = 1. This 
corresponds to a uniform flow throughout, and, while 
it is a valid solution, it is not of particular interest 


here and, hence, is disregarded. From Eqs. (6b), 
(Sa), and (9a), it follows that 
T ¥-1 u ” 
a lL -+ 5 Me (1 — V?) 
Hence, Eq. (12a) can be written in the form: 
— |] = n . ; 
f 1+ 21 — Vv) | Val 
t= 4 _ arn ene + C, 
«Mo V — 1)(V — a) 
(12b) 


For general values of » this integral must be evaluated 
numerically to give & = &(V). 

To complete the description of the flow given by this 
it is convenient to determine the other flow 
Thus, 


solution, 
variables as explicit functions of V. from Eqs. 


(2) and (9a), 
p/po = Uo/u = 1/V (13a) 
As before, 


7 l 
ME te a. ok 


- —— Mo%(1 — V2) (13b) 


*« = 1.36 for air with y = 1.40 
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Hence, from the gas law — [Eq. (le)], 


a Wat 7s 
po V 


Furthermore, it will be of interest to determine the en- 
tropy S and Mach Number J/ as functions of V. In 
general, for a perfect gas, 


(13c) 


MMe. — v9 | 


S — So = c In (T/T) — R In (p/po) 
Substituting Eqs. (13b) and (13c) in the above, one 
finds 
S — So 


= te {hi 42 
Co 


If the sonic velocity be denoted by a, then 
Mu “ (“*) =, (F) 
My wa) To 


M y¥-1 mae : 
—=V|1 — Mo(1 — V? (13e) 
Mo | > ] 


— 1 
Mer(1 — v9) | v7-*ha3a) 


or 


The solution given by Eq. (12b) is essentially the one 
obtained in a restricted manner by Becker* for n = 0 
and by Thomas‘ for n = '/s. For completeness the 
solutions for nm = 0 and n = 0.768, the experimentally 
determined value for air, will be given explicitly. 


Constant Coefficients of Viscosity and Conductivity 


If the coefficients of viscosity and conductivity are 


assumed constant, then x = 0 and Eq. (12b) becomes 


Bact 
(V — a)* 


where C,’ is an arbitrary constant. 
V = last— —o, itis necessary from Eq. (14a) that 
a <1. This, in turn, implies that 1/4) > 1—that is, 
that the flow, in order to be other than the trivial com- 
pletely uniform must be initially supersonic. 
Thus, in this formulation the necessary condition of an 
initial supersonic Mach Number for the existence of a 
shock wave is mathematically obtained without the 
use of the second law of thermodynamics. 

For C;’ negative, Eq. (14a) yields a monotonically 
with V becoming infinite and 7 
becoming negatively infinite as——> +. Such a flow, 
to the authors’ 


= C,’ exp [x(1 — a) Mok] (14a) 


Since by definition 


flow, 


increasing velocity, 


the mathematical existence of which, 
knowledge, has not been previously established, would 
be completely supersonic and would correspond to an 
expansion wave. However, it evidently can be valid 
only in a finite region wherein the temperature is posi- 
Moreover, it will be shown that such a flow would 
there- 


tive. 
involve a continuous decrease in entropy and, 
fore, would not be obtained physically even in a finite 
It is interesting to note that such a solution 
at plus infinity, 
Thus, 


region. 
would 
where both V 


nonuniform conditions 


and d/dé would become infinite. 


give 








even when C, = O and Eq. (6b) is valid, mathemati- 
cally possible solutions with nonuniform conditions at 
plus infinity can be obtained. 

If C;’ is determined such that the origin of the coor- 
dinate system is located where the velocity profile has 
an inflection point—that is £ = 0 where d?V/d& = 0 
then C,’ is positive and Eq. (14a) becomes 


1 — Va 
(Va — ee 


or, in another form, 


. l , |(~ = *) i — 7: | (140) 
¢= —— n : e ) 
(1 — a)Mp V-—a /1 -—Va 


From Eq. (l4c) and Eqs. (13a)-(18e) the entire 
characteristics of the flow can be determined. Exami- 
nation of the variation of V, p, p, 7, S, and M with é 
indicates that the flow described by this solution is the 
well-known compression or shock wave, which satis- 
fies the Rankine-Hugoniot values at plus and minus 
infinity. Figs. 1 to 3 show the distribution of V, p, 
and S with ¢ for Mp = 1.1 and for Mo = 2. 

The entropy distribution given by this solution may 
warrant a brief discussion. It will be observed from 
Fig. 3 that the entropy increases from zero at minus in- 
finity to a maximum at the center of the wave (£ = 0) 


1— 


V 
(V a) exp [k(1 — a@) Mot] (14b) 
tc 


and then decreases to a positive value at plus infinity 


D0 == = 7 
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This maximum in the entropy at the center of the wave 
can be demonstrated analytically as follows. In gen- 


eral, 


(15a) 


a Pe 
dx y—1Tdx udx 


From Eqs. (9a), (10), and (13b) with m = 0 and C; = 
0, one finds that Eq. (15a) may be written as 


2 lL es 





i 4 - [(¥ = 1h(V - a)(V? — @) 
y+ 1 «RM, dé - Vy? T/T» 


(15b) 


Eq. (15b) indicates that the entropy has stationary 
values at three values of V—viz., V = 1, Va, a. Itis 
clear that the stationary value of S at the center of the 
wave, where V = Va [cf. Eq. (14b)], is a maximum, 
since Eq. (15b) shows that dS/dé is negative fora < V 
< Va—i.e., for V downstream of the center. It may 
at first sight appear that the recovery of mechanical 
energy on the downstream side of the center of the wave 
thus indicated by this solution would violate the second 
law of thermodynamics and thereby invalidate the 
solution. However, the second law applies to an entire 
system—that is, to the end points of the flow—and 
permits energy recovery in separate sections thereof. 
The negative entropy gradient here might also be in- 
terpreted as indicating physical effects that are not 
taken into account by the governing Eqs. (la)—(le). 
For example, the gas molecules in some regions may not 
be in thermal equilibrium.’ 

It is clear from Eq. (15b) that the expansion-wave 
solution discussed above must be discarded on the basis 
of entropy considerations. For this solution the veloc- 
ity increases monotonically with §—that is, V > 1. 
It therefore follows from Eq. (15b) that dS/d& < 0 
everywhere, which is physically impossible. 

The solutions given by Eqs. (14a)—(14c) are valid 
solutions of the integral in Eq. (12b) provided that a 
<landn=0. Itisof interest to consider the case of 
n= 0 but a = 1—that is, Jo = 1. In that case Eq. 
(12b) gives 

t= {tn Cy"(V — 1) + ws | (16) 

' K _ 1— V 
where C,” is an arbitrary constant. In order that V = 
| at minus infinity, it is necessary that C,” > 0, for, if 
‘" < 0, then it would be necessary that V < 1 for & 
real, and it would then follow from Eq. (16) that as 
V—>1,&—- +o instead of > —. For a positive 
value of Cy”, Eq. (16) represents a flow in which the 
velocity increases monotonically as £ increases. This 
is evidently a special case of the expansion-wave flow 
discussed above and must be discarded on the same 
grounds (entropy) as that solution. Thus, the only 
physically possible solution for the case Mp = 1 is the 
solution V = 1 everywhere. 
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It is thus shown that a flow that is sonic at minus 
infinity will remain sonic without change. More gen- 
erally, it has been proved here that, if a one-dimensional 
flow is subsonic or sonic at minus infinity, then the flow 
will be uniform throughout the entire range of x. This 
is true for the cases of both uniform and nonuniform 
conditions at plus infinity, since in the latter case, the 
additional term (proportional to C.) in Eq. (6a) van- 
ishes at minus infinity, so that the mathematical argu- 
ment upon which the above statement is based remains 
valid. 


Variable Coefficients of Viscosity and Conductivity 


For n + 0, the integral in Eq. (12b) can be evaluated 
numerically for any given m and Mo. The arbitrary 
constant C, may be determined, as in the case n = 0, 
by placing the origin (£ = 0) at the inflection point 
(d?V /d& = 0) in the velocity-distribution curve. From 
this numerical solution for — = &(V), it is possible to 
determine the variation of the other flow variables by 
using Eqs. (13a)—(13c). 

The experimental value of n for air is 0.768. The 
numerical evaluation of the integral in Eq. (12b) has 


been carried out for this value of and for Mp = 1.1 
and 2.0. The results are shown in Figs. 1-3, along 
with the results for » = 0 for ready comparison. As 


might perhaps have been expected, the increase in vis- 
cosity and conductivity through the shock wave de- 
creases somewhat the absolute value of the gradients 
of the flow variables in these two numerical examples. 


SHOCK WAVE THICKNESS 


Mathematically, the changes in the flow variables 
through a compression shock occur over the infinite 
distance from minus to plus infinity. Inspection of 
Figs. 1-3 will indicate, however, that in reality the 
major portion of the changes occur in a small region at 
the center of the wave. This small region is called the 
shock-wave thickness. The question which must be 
considered is whether the shock-wave thickness is of 
the order of magnitude of the mean free path of the gas 
molecules. If it is, then the ability of the above equa- 
tions, which are based on the continuum hypothesis, 
to describe the structure of the shock wave is doubtful. 
Previous authors?~‘ have indicated that only for weak 
shocks (Jp < 1.3) is the continuum theory valid. 
However, because of the rather arbitrary manner in 
which shock-wave thickness has been defined, it was 
thought advisable to compare the results given by the 
above analysis according to the various definitions. 

Nondimensionalizing the space variable x by means 
of the mean free path (A)* of the gas molecules at 
minus infinity is convenient for the consideration of 
shock-wave thickness [cf. Eq. (9a)]. The thickness is 
then given as a pure number—namely, the number of 


* For air with 7) = 519°R., po = 0.002378 slugs per cu.ft.; 
that is, for standard sea-level conditions, 4) = 3.46 K 10~%in. 
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mean free paths at minus infinity. If this number is 
large in comparison with unity, then the continuum 
hypothesis is valid; if this number is of the order of 
unity, then the continuum hypothesis is no longer a 
sufficiently accurate basis of determining the flow char- 
acteristics. For purposes of a more refined appraisal, 
the thickness may be calculated as a multiple of the 
local mean free path of the gas molecules at the shock 
wave rather than of the mean free path at minus infin- 
ity. 

The point of view that may be taken here regarding 
the mean free path is that this quantity can be deter- 
mined experimentally from the measurable quantities 
viscosity, density, and temperature. Thus, A, in gen 


eral, may be consistently defined by the equation 


(cf. Eq. (9b) } 
A = (1/6) (uV r/pV SRT) (17a) 
In particular, Ao will be 
No (1/5) (uo r/po SRI l7b) 


This point of view is consistent with Eqs. (9b) and (9c 
regardless of the value of 7.* 

The various definitions of shock wave thickness will 
now be discussed and compared. Prandtl'® has given 
a definition of shock-wave thickness ¢,; shown in Fig 


fa. It may be stated mathematically as 


Noth = —(uo — Uy) /(du/dx) mar (1Sa 


Since the origin of the space variables x and & has been 
placed at the inflection point in the velocity distribution 
and since V 
of V 


= aat plus infinity, this becomes, in terms 
and &, 


ty -(] a)/(dV/dé)_ (1Sb) 

Caylor and Maccoll’ have given another definition of 
shock-wave thickness f-—namely, the distance in which 
a certain portion P of the total velocity change through 
the wave takes place (see Fig. 4b). They selected P = 
0.80 and, according to this definition, gave a simple 
approximate formula for the thickness of a weak wave. 


Constant Coefficients of Viscosity and Conductivity 


The shock-wave thickness ¢,’ for » = 0 according to 
Prandtl’s definition will be calculated first. Eq. (14b) 
gives V = Va for ¢ = 0, and thus, from Eq. (10) with 


n = Oand C, = 0, Eq. (18b) becomes 


hh’ = (1 — a)/«eM,(1 — V/ a)? (19) 

* If the molecules of a gas are assumed to be spheres of diame- 
ter o and mass m, then the mean free path can, by kinetic theory, 
be shown to be A = m/+/2mrpo*?. Hence, if o were constant, then 
\ would be a function only of the density p at any point, and ac- 
cording to Eq. (17a), « would be proportional to 7'/? implying 
1/2, Thus, by choosing n ¥ '/2, one is implicitly takinginto 
account the variation of molecular diameter with temperature. 
In fact, with the equations used here, ¢ ~ T “~- 2/4, More- 


over, if nm # 1/2, then X is implicitly considered to be a function of 


1 = 


T as wellas of p. 


AL 
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Values of ¢;' for a range of values of J) are shown in 
Fig. 5. 
The thickness f,’ according to the Taylor-Maccoll 
definition for » = 0 is calculated for any value of P as 


follows: Eq. (14c) gives (ef. Fig. 4b) 


te’ & + & = 
l In l ate - a\ (20a 
K(1 — a)Mo l= VV -— a/ 
where 
i+ P | P 
V/V; = a = (l — a) = (] Ta =U a 
i—P r 
Vo act+ (l— a) = (l+a)— (1 — a) 
9 9) 2 


or, finally, 
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i ae 
be In (20b) 
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Values of t,’ have been determined for P = 0.80, 0.90, 
0.95 for a range of Mp» and are shown in Fig. 6. 

Examination of the values of ¢,;’ and ft’ shown in 
Figs. 5 and 6, respectively, indicate that the Taylor- 
Maccoll definition imposes a somewhat less severe re- 
striction of the continuum hypothesis than does the 
Prandtl definition. However, if one rather arbitrarily, 
but reasonably, requires that ¢ > 7 for the above govern- 
ing equations to describe a normal shock wave ade 
quately, then according to any definition it would be 
required that My < 1.3. Thus, only for weak waves 
would these equations describe the actual conditions 
ina shock wave. 

Two refinements in the calculation of these thick 
nesses may be made: consideration of the variability 
of the coefficients of viscosity and conductivity, and 
determination of the thickness in terms of the local 
mean free path at the wave. 


Variable Coefficients of Viscosity and Conductivity 


For determining shock wave thickness ¢” when the 
coefficients of viscosity and conductivity of the gas are 
considered variable, only the Prandtl definition will be 
used, since it is easier to apply here than the Taylor- 
Maccoll definition and since, as has been seen above, 
the two definitions lead to similar results. 

In order to apply the Prandtl definition in this case, 
it is necessary to calculate the value of V at the origin 


of the axis. By differentiating Eq. (10) (C. = 0) with 
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respect to £, setting d?V/d& = 0, and eliminating dI’/- 
dé, one obtains the following equation for the value of 
|’. -that is, for the value of V at & = 0, 


—. 
(» - ‘) VA—n (1+ a)V3 4 


: ayt+] 
( : tn)ave- 2? =0 (21) 
y¥-—I!1 27-1 


The value of V. may be determined from this quartic 
equation for any given » and a(or J/o) by Newton’s"' 
(approximate) or Ferrari's (exact) method.'’** From 
this value and from Eq. (10) (C. = 0), (dI’/dé); =o 
may be determined, and then /,” can be calculated from 
Eq. (18b). 


The calculations outlined above were carried out for 
n= '/, and n = 0.768 for a range of AJ) values. The 
results are shown in Fig. 5. As might be expected, 
the wave thickens as » increases. 
pointed out by Thomas,’ for m = '/s, the thickness 
However, 


Moreover, as was 


reaches a minimum finite value as JJ) > ©. 
for n > '/s, the thickness decreases to a minimum 
(about 1.5 mean free paths at Jp = 5 when = 0.768) 
and then increases without limit as IJ) — ©. This 
can be demonstrated analytically in the following man- 
ner: Eg. (10) (C2 = 0) can be written, with the aid of 
Eq. (13b), as 
dV K(V — 1)(V — a) 


— : (22a) 


ae” fi 7-1 2%, Joe | 
a4 1— v)| Mo” 
lana a , 


From Eq. (12a), it is clear that as My > ©, a — 
(y — 1)/(y + 1). Moreover, from Eq. (21), V. ap- 
proaches a finite value, depending on m, as Mp > @. 
Thus, 
j dV K(V. — 1)(Ve — a@) 
lim ~) = or | 
lo ES 7E =0 , 2 ai | 
Vv t V. k a (h- v2) | on 
“ ( (22D) 
ifn = '/s | 
=-=@Q f2> "/s; | 
—->o, ifn < '/s 


’ 


It follows then from Eq. (1S8b) that, as JJyp — ©, 
the shock-wave -thickness given by this analysis de- 
pends strongly on the value of m; it will have a finite 
value if m = '/s, will approach zero if m < 1/2, and will 
approach infinity if m > ‘/2. Thus it is seen that, if 
n > '/, (as is actually the case for air) it may not be 
quite true that the stronger a normal shock wave, the 
thinner it will be. However, the thickness increases 
slowly for n = 0.768 as Mo increases beyond the value 
(My) = 5) at which the thickness is a minimum, while 


phenomena such as dissociation, ionization, and con- 


* One must, of course, find the (real) root that is between V = 


land V=a 
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densation not considered here may then become signifi- 
cant. 

Examination of Fig. 5 indicates, therefore, that con- 
sideration of the variability of the coefficients of vis- 
cosity and conductivity does not appreciably change 
the Mach Number at which the continuum hypothesis 
must be considered of doubtful validity. If » = 
0.768 and ft,” is required to be greater than 7, then Mp 
must be less than 1.3. 


Thickness in Terms of Local Mean Free Paths 


If one is interested in the actual dimensions (say in 
inches) of the shock wave as given by the two defini- 
tions, the previous results are adequate. Thus, for 
particular conditions at minus infinity (1p, po, and 75), 
the thickness of the wave can be determined by calcu- 
lating Ao from Eq. (17b) and by multiplying the ap- 
propriate value of ¢ from Figs. 5 or 6 by this value of 
Ao. However, to examine the validity of the contin- 
uum hypothesis, a more suitable appratsal of the shock- 
wave thickness may be obtained by determining this 
thickness as a multiple of the (local) mean free paths 
in the wave. 

A good approximation to this local value may be 
found by taking the arithmetic mean of the mean free 
paths of the gas molecules at minus and plus infinity. 
Thus, 


ty — 2t [1 + (Ay Xo) | (23a) 


where /, is the thickness of the wave in multiples of the 
mean free path in the wave, ¢ is the thickness in multi- 
ples of the mean free path at minus infinity, and \, is 
the mean free path of the gas molecules at plus infin- 
ity. From Eqs. (9c), (17a), and (17b), it can be shown 
that 


(1/2) | (23b) 


Ay Ao = VT; 7" 


Observing that V; = a [cf. Eq. (12b)], Eqs. (13b) and 
(23b) may be used to write Eq. (23a) in the form 


2t 
h=— -: 
» y¥— l : [m — (1/2)] 
1+a|]1+ —— M1 — a’) 


= K(M, n)t 


where K (Mo, 2) is a multiplicative parameter depend- 
ent only on Mo and m. Obviously, 4, as defined by 
Eq. (23a) will, in general, have a value between ¢ and 
2t. Applying this parameter to the thickness ¢,” for m = 
0.768 results in the thickness shown in Fig. 5. 


From a consideration of Fig. 5 the conclusion may be 
drawn that for ¢ > 7, Mo must be less than approxi- 
mately 1.3. Thus, considering the variation of mean 
free path through the wave does not significantly af- 
fect the conclusions based on the analysis in terms of 
the mean free paths at minus infinity. 
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TERMINATING SOLUTIONS 


If C, * 0, then it can be seen from Eq. (6a) that, at 
x —> ©, at least either u or 7 becomes indefinitely large 
Thus, the region of the flow must in this case be finite 
in the direction of the positive x axis. In other words, 
the flow must be terminated at some point (denoted by 
£ in the & plane), at which certain conditions deter- 
mining C2 are prescribed. If, on the other hand, physi- 
cally significant conditions are to be prescribed at 
plus infinity, then these must be the Rankine-Hugoniot 
conditions, and C, will be zero. This implies that a 
necessary and sufficient condition for uniform flow be 
hind a shock wave is that the Rankine-Hugoniot rela- 
tions between the thermodynamic variables before 
(—o) and behind (+ ©) the wave be satisfied. 

For the purpose of illustrating the nature of solutions 
of the flow equations for which C, ¥ 0, the coefficients 
uw and k will henceforth be assumed constant, implying 


n= 0. 
Let 
z = exp {2KyMoEt/(y + 1)} (24a) 
Then Eq. (10) can be written in the form: 


1V 1 ein Pe 
Vs‘ ee 1\(V — a) - : - =0 
dz 27 Y Uo- 


(24b 


A complete numerical solution of Eq. (24b) can be ob- 
tained by prescribing the Mach Number Mp at £ — 
—o and by prescribing the values of two variables’ 
(say, 71/T>) = 71 and %/u% = V;) at some finite point 
to be chosen as the origin £ = 0. In the 2 system, the 
range — © < £ < Ocorresponds to the range0 <2 <1. 
The values of the dimensionless constant C2/o” can be 
obtained in terms of Mo, 7:, and V, from Eqs. (6a) and 
(Sa). One thus finds 


c. l l - 
= V;' -(r,; — 1) (24e 
Uo~ $ 7 > 1M? 


Three types of flows can be mathematically ob- 
tained, depending on whether C2/1o is zero, positive, or 
negative. If C./uo? = O, then the quantity // 
(u2/2) + cpT, which in adiabatic flow represents the 
enthalpy at the stagnation point, will be constant. 
Moreover, it follows from Eq. (24c) that, if the Ran- 
kine-Hugoniot relations are satisfied at point ‘1’’ for 
two of the thermodynamic variables (e.g., for V aud 
for T/T»), then C2/uo? = 0. If Co/uo” > O, then // in- 
creases in the direction of flow, while if C,/mo? < 0, then 
H decreases. It is interesting to note, in fact, that in 
the one-dimensional theory an increasing H downstream 
is mathematically consistent with the law of the conser- 
vation of energy as stated by Eq. (Ic), as well as with 


* If p/po and p/p» are given, then T/T» and u/uo can be de- 
termined by means of the continuity equation and the equation 


of state. 
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the second law of thermodynamics (cf. discussion of 
results below). 

With the numerical value of C./m9? determined, the 
first-order differential Eq. (24b) can be solved numeri- 
cally, without much difficulty, by standard methods,* 
starting from the point (¢ = 1, V = V;) and proceeding 
backwards to z = 0 and forward to some finite value 
ofz > 1. (Asacheck it should automatically be found 
that V = latz = 0.) 

The case C:/m? = 0, leading to compression shocks 
satisfying the Rankine-Hugoniot relations before and 
behind a front, has already been treated here in detail. 


*Cf., for example, reference 11. Also cf. method of Taylor 
series described below. 
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To indicate the nature of solutions in the case C2/u? + 
0, two numerical examples will be discussed here: one 
in which C:/u2 > O and one in which C3/u? < 0. 
In both cases the values V; and 7; will be chosen so that 
the second law of thermodynamics will not be violated 
at least with respect to the points § ~ — and — = 0 
i.e., so that S; > So.F 


Numerical Examples 


») 


Example I.—Suppose p~:1/po = 4 and 7,/T) = 1 = 2. 
Also suppose Mp = 2.{ Then from Eqs. (le) and 


(2), it follows that 2;/m% = Vi = '/2. Eq. (24c), with 
y = 14, then gives 
C2/uo? = 1/4 


From Eq. (6a), it will be found that this implies that 1 
will increase, at £ = 0, to 1.222 times its original value 
atmé— —o, 

Eq. (24b) was solved numerically by starting from 
the given values of V and its derivatives at s = 1 and 
expanding about this point in a Taylor series. This 
expansion was carried out for the limited region 0.5 < 
z < 1, so that only a few terms in the series were needed 
for sufficiently rapid convergence. For the next region, 
0.2 < s < 0.5, another Taylor series was used, this time 
about the point z = 0.5. This procedure was continued 
toz.= 0. Near z = 0, however, it was necessary to 
use much smaller intervals of z in order to obtain satis- 
The expansion was similarly 


factory convergence. 
The final 


carried out forward for — > 0, or z > 1. 
results are shown in Figs. 7 and 8. 
Example II.—Suppose V; = 0.6, 7 = 1.5, and My = 
2 (y = 14). 
Proceeding as above, one finds 
C2/uo? = —0.0075 


indicating a decrease in //, at — = 0, to 0.993 times its 
value at & Eq. (24b) was solved numerically 
by the Taylor series method described above, and the 


>—o, 


results are shown in Figs. 9 and 10. 

The entropy (S — So)/c, was calculated for both 
examples, and it was found that for the region investi- 
gated the quantity S — So was everywhere positive, 
so that the second law of thermodynamics was not 
violated. 

It can be seen from Figs. 8 and 10 that the solutions 
with C, ~ 0 yield compression waves followed by large 
expansions.** These compressions and expansions all 


T This condition is imposed so that the boundary conditions 
correspond to physically realizable flows. 

t It may be objected that, as seen previously, the case Mp = 2 
may give unreliable results when treated by the continuum the- 
ory. The purpose of the present examples, however, is only to 
indicate the nature of the solutions obtainable, so that the par- 
ticular value of Mp» chosen does not matter. 

** Some rough calculations have indicated that the pressures 
in Figs. 8 and 10 must continue to drop considerably before any 
compression can occur. Moreover, it can be shown that dp/dt — 


—wasti—- om, 
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occur in a small region—namely, a region of the order 
of magnitude of several mean free paths—although for 
a smaller initial Mach Number the extent of this region 
increases. The physical significance, if any, of these 
types of solutions remains open to question. 


CONCLUSIONS 


Based on the one-dimensional analysis of a continu- 
ous, viscous, heat-conducting fluid given here, the fol- 
lowing conclusions can be stated: 

(1) The Eq. (u?/2) + ¢,7 = constant is valid, only 
if the Prandtl Number has a constant value of °/4. 
Moreover, this equation is a complete integral of the 
energy differential equation only under the condition 
of uniform flow at plus, as well as at minus, infinity. 
The particular value P, = ; here is based on the 
Navier-Stokes equations and is therefore subject to the 
physical limitations of these equations. 

(2) Mathematical solutions representing monotonic 
expansion waves can be obtained. These solutions 
violate, however, the second law of thermodynamics, 
even in the limited region in which the temperature 
would be positive. 

(3) Under conditions of uniform flow at =~ and for 
a Prandtl Number of */4, simple closed-form solutions 
of the flow equations, which represent monotonic com- 
pression shocks satisfying the Rankine-Hugoniot con- 
The thick- 


ness of such waves, which, in general, are of the order 


ditions at the end points, can be obtained. 


of magnitude of several mean free paths, will decrease 
as the Mach Number .\/) increases, provided that n < 
'/5. If, however, m > '/»s (as for air), then the thickness 
will diminish with \/ until it reaches a minimum (about 
two local mean free paths when J/) = 5 if 2 = 0.768), 
beyond which the thickness slowly increases to infinity 
as Mp increases to infinity. The continuum hypothesis 
probably gives inaccurate results for those cases in 
which the thickness is small in comparison with the 
mean free path (cf. Fig. 3.) 

(4) The Rankine-Hugoniot conditions are necessary 
and sufficient conditions for the flow to be uniform at 
both plus and minus infinity, 

(5) A flow that is either subsonic or sonic at minus 
infinity will remain uniform throughout the entire 
field of flow. 

(6) Mathematical solutions of the flow equations 
can be obtained for which the quantity (u?/2) + c¢,T 
varies, 
infinite regions, since the flow variables become infinite 
as § > +o-. 


waves followed by expansion waves in small regions of 


Such solutions are applicable only to semi- 
These solutions indicate compression 
the order of several mean free paths. The physical 
significance, if any, of such solutions remains, however, 
open to question. 


(7) Phenomena observed in actual steady normal 
shock waves which cannot be explained by any of the 
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types of solutions given here must be due to either multt- 
dimensional effects or to effects not contained in the 
hydrodynamic continuum hypothesis of a gas with 
constant coefficients of specific heat and following the 


ideal gas law. 
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On Shock Waves in Inhomogeneous Flow 


MAC C. ADAMS* 
Cornell Unwersity 


SUMMARY 


In connection with certain problems involving mixed subsonic- 
supersonic flow fields, and especially problems of normal shock 
waves in channels, it may become important to determine the 
configuration of a nominally ‘‘normal’’ shock wave in a slightly 
inhomogeneous steady field of flow. In this note it will be shown 
that the position assumed by the shock wave, its strength, and 
the details of the flow behind it can be calculated by a small per- 
turbation theory, assuming that the velocity disturbances in the 
supersonic region are small compared to a parallel stream velocity. 
This calculation will be carried out in detail for the two-dimen- 
sional case—i.e., the case of plane disturbances of a parallel 
stream. The extension to axisymmetric cases would seem to 
offer no essential difficulties 

It will be necessary to account for first-order vorticity in the 
subsonic region; this can be done by means of the differential 
equation derived by Sears.' Detailed results will be presented 
for two different typical disturbance profiles. 


INTRODUCTION 


: bes ANALYSIS BEGINS with the assumption that the 
undisturbed flow consists of an unbounded parallel 
stream traversed by a normal shock wave. This, of 
course, represents a simplification, for a diverging 
channel would be necessary to produce a normal shock. 
It is assumed here that the divergence is negligibly 
small and that the walls are so far away from the region 
under consideration that they have negligible influ- 
ence. 

The treatment of a bounded region of flow would be 
similar to that presented here, with images or similar 
mathematical devices introduced to satisfy the bound- 
ary conditions at the walls. 


ANALYSIS 


Let the axis x = 0 lie along the normal shock wave of 
the undisturbed flow so that it divides the flow field 
into two unbounded regions—supersonic in the region 
x < 0 and subsonic for x > 0, with Mach numbers 
Mh > land Mz < 1, respectively. Thus, 1/2 is a known 
function of 4, and y the ratio of specific heats. The 
disturbed shock wave is assumed to deviate only 
slightly from the y-axis, and conditions appropriate to 
the upstream and downstream sides of the disturbed 
(See Fig. 1.) 


wave will be applied atx —~ +0. 
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tical Engineering. 


It will be convenient to discuss separately (1) the 
conditions in the supersonic and subsonic regions of 
flow and (2) the passage through the shock wave. 


Supersonic Region (x < 0) 


Let 1%, 7; denote the disturbance velocity components 


in this region, and assume that they are small compared 


to the stream speed U;. They must then be related by 


the well-known equations of linearized supersonic 
flow—i.e., 
u(x, y) = F’(x — my) | 
= (1) 
n(x, vy) = —mF'(x — my) 
where m denotes V \/,;* — 1 and the prime represents 


differentiation. For example, if the disturbance is 
produced by a solid boundary defined by the curve y = 
h(x), (« < 0), then the linearized boundary condition is 
v(x, h) = U,h'(x) and the streamwise velocity compo- 


nent is given by 
u(x, y) = —(1/m)Uih'(x — my) (2) 


Crossing the Shock Wave 


The disturbance velocity components just upstream 
and downstream of the shock wave are related by the 
Rankine-Hugoniot shock-wave equations. It will be 
shown, however, that, under the assumption of small dis- 
turbances, the streamwise velocity component u im- 
mediately downstream of the shock wave is uniquely 
determined, in the first order, by ™. 

Consider any local point along the disturbed shock 
(see Fig. 2). From the figure, one deduces the following 


approximate relationships: 
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Fic. 2. Diagram showing angles at shock wave. 


OB = U1 + 1, AC = (U2 + U2)0 
1/c = AB/AC = (OB — OA)/AC 
6 = AC/OA 
Thus, 
6/0 = [(Ur + m)/(U2 + u2)] — 1 (3) 
Another relation for 0/o can be obtained from the 
shock-wave equations :” 
y+ 1 MT? 
- 


M? cos? vam ] 


tan o 


= aaa (4) 
tan 0 

Here M denotes the local Mach Number upstream of 

the shock wave. Let /’ denote WM — 14; then M’ << 


1 and, to first order, Eq. (4) becomes 


6 ae Me = 4 (y+ z Toe - 
7 — 7 a vo 
1+ —— My; (1 + —* us) 
Moreover, in the same approximation, 
ea Se el (6) 
a4+a ay ay a," 


The quantity a’, representing the disturbed velocity of 
sound, is found by use of the energy equation and is 
given approximately by 


a’ — —I[(y — 1) 2) Many (7) 
Combining Eqs. (3), (5), (6), and (7) yields 
U2 = —(U2/U;)u; (8) 


It is observed that this result is the same as would be 
obtained for a normal shock with only a streamwise 
velocity disturbance. 
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Subsonic Region (x > 0 


It now remains to relate the vertical velocity per. 
turbation in this region to w2, and therefore %. Then 
the shock angle can be calculated at any point along 
the axis x = 0, and the shock position easily follows by 
integration. 

As will be shown later, the vorticity in the subsonic 
region is of first-order magnitude, and, consequently, 
potential flow will mot exist. The stream function 
y is therefore introduced and is defined by 


vy = (p/p2)u, —z = (p/p2)v (9 


where u = U,; + u’ andv =v’. The quantity p repre. 
sents the density, and the reference density p: is con- 
veniently chosen to denote the density that the free 
stream would attain in passing through a normal shock. 
From here on, all symbols without subscript will denote 
conditions downstream of the shock, and the subscript 
2 will identify conditions of the free stream after pass- 
ing through a normal shock—1.e., conditions in the 
subsonic part of the undisturbed flow. 

If the velocity perturbations are considered small as 
compared to U; and if the stream function is written 


(x, vy) = Uoxy + v(x, 9) (10) 


then the following differential equation is obtained:! 
By.’ +,’ = -[1 + (vy — 1M, (11 
where 82 = 1 — M,.? and Q = 2,’ 


shown, in reference 1, by consideration of the energy 
equation and equation of state of a perfect gas, that 


vy’ a B2u’ a [(S _ Se) R] Us, —y,' —_ vy’ 


The quantity S denotes entropy, and R is the gas con- 
stant. 

The appropriate solution to Eq. (11) can be con- 
structed by making use of the following transformation, 
which puts Eq. (11) into the form of Poisson’s equa- 
tion: 


— u,’. It is also 


X =x/8B, Y=y 
y'(x, vy) = f(X, Y) 


The solution consists of the particular integral involv- 
ing 2 and a complementary function that may be con- 
sidered as representing a distribution of sources along 
the Y-axis: 





_ Ore re .. = 
SX, Vy = — 32 ff "0G apinlcx — 98 + (¥ — nas da + 5 fet) tant a (12) 
Ar —_ Qn J_.« X 


where 


Co = ] + (y as 1) M,?. 


Differentiating and transforming back to the real plane yields 


S— So - C, = o Q é, o 
pul Ste ~ 3 fae [” ee ant ef, —— 
R 2r Jo ~o (x/B — £)? + (y — »n)? 22 J — a (x/B)* + (y — 9)? 


a i -/B)g( 
F (x/8)g(n) (13) 
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ON SHOCK WAVES IN 


The vorticity 2 and the entropy S must now be deter- 
mined. This will yield relations for the perturbation 
velocities immediately downstream of the shock, since 


u’(0, vy) = ue, v'(0, y) = v% (15) 


The source distribution function g(n) will then be read- 
ily determined, since the second integral of Eq. (13), for 
x = 0, becomes® 

(1/2)g(y) 
subject to the restrictions that the function be con- 


tinuous and bounded. 
To determine the vorticity, use is made of Bjerknes’ 


theorem: 

q X 2 =T grad S (16) 
where g denotes the vector (U2, + u’) + jv’ and T 
represents the temperature. Eq. (16) can be written, 
in part, 

(Uz + u')Q = T(OS/dy) (17) 

and to first order 

Q = (T2/U2)(0S/Oy) (18) 


Now, since the inclination of the flow, v’/ Us, is assumed 
extremely small, the terivative 0S/Oy constitutes an 


approximation to the gradient: 
oS/Oy cdl grad S | 


grad S | is approxi- 
It can therefore 


and in small perturbation flow | 
mately constant along stream lines. 
be seen that the vorticity can be considered here to be 
nearly constant along any stream line. 

Since the entropy can be determined only immedi- 
ately downstream of the shock, a further simplification 
must be made in order that the vorticity can be calcu- 
lated. This is to consider Q(x, y) = Q(y). This ap- 
proximation is consistent with oe. first-order theory 
being employed here, because the inclination of any 
stream line relative to the x direction is small, as men- 
tioned above, and the error introduced in the calcula- 
tion of velocity components at x = 0 is of second order. 
This is analogous to the approximation made in the 
familiar Prandtl wing theory regarding the location of 
the trailing vortex sheet behind the airfoil; it is well 
known that this causes only higher-order errors in the 
induced velocities at the wing. 

The entropy at any point along the shock, on the 
downstream side, can be found by use of the equation® 


M? cos? « — 


+1 
¥- 


‘Ae — 1)M* cos? a + =| (19 
: : 9) 
y+ 1 (y + 1)M? cos? o , 


INHOMOGENEOUS 


Co Q(E,n) (v/B — &) 
y! = d 
of me B-b?+ 0-9" 2 = B 


FLOW 687 


= g(n)(y — 1) dn 


rr (14) 
(x/B)? + (y — 9)? 
where 5S; is the entropy of the free stream. Expanding 


Eq. (19) in a Taylor series about 17 = M, and o = 0, 


one has, approximately, 


S — S pe 27/( M\? — i uy (20) 
R 2yM -— (vy - 1) U; 
and therefore 
ay) = i A 2 foe (21) 


[2yMi? — (vy — 1)] M2? U; oy 

In this paper calculations for the subsonic velocity 
perturbations will be made only immediately down- 
stream of the shock; however, Eqs. (13) and (14) could 
be used to determine the velocity at points far down- 
stream. 


APPLICATION TO SINUSOIDAL DISTURBANCE 


As an example of the use of the equations developed 
above, the configuration assumed by a ‘“‘normal’’ shock 
wave in a field of small sinusoidal disturbances will be 
determined. Such a disturbance field might be gen- 
erated, for example, by a wavy wall whose shape is 
described by 


h(x) 


= Ecosmx, —~» <x < — myo 
In this case 


—(xE/m)U, sin rmy)\ 99) 


u,(—0, y) = 
(rE/m)U2 sin rmy § - 


u2(+0, y) 

The part of the shock wave to be investigated here is 
that part far up along the y-axis. In other words, the 
is located so that the dimension 


origin of coordinates 
—yo is extremely large, and Eqs. (22) can be assumed 
to apply for all values of y. 


From Eqs. (20) and (21), one can write immediately 
| z 


$~s (M,2 — 1)" + | 
— = rit sin mmy 
R 2yM,? — (y — 1) “ } (23) 
= —() sin mmy, say 
2(M\? — 1)?U2r*E 
Q(n) = -= ot wf cos mmn (24) 


[2yM,? — (vy — 1)] M.? 
= —C, cos mmm, say 
Setting x = 0 in Eq. (13) and carrying out the inte- 
gration gives the source distribution 


E OoCi\ . : 
g(y) = (2n6°U: + 20,U2 - 2 sin rmy (25) 
m n 


Then, by use of = (14) and after simplification, 


m — xBMY{ M2 — [Cy + 3)/2) 
U. as 1) /2) My? 


= ¢ cos amy (26) 
okt m 
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M, \w 
Fic. 3. Ratio of shock wave to wall amplitudes for sinusoidal | 
flow disturbance, as function of stream Mach Number, for air 
(The dashed curve shows the results of a calculation in which \ Ni WI 
vorticity is neglected.) \ r ua 
laa al 
For future use the quantity in brackets will be denoted ! 
by the symbol A. When 1/, = V(y + 3)/2,A =0 2 0 | 2 - 
and the flow is horizontal immediately downstream of | ™N Exe 
| oe ; | me 
the shock. It will be seen later that this same Mach € | 
Number produces other interesting phenomena. | | 
= rs + 
Now that all the velocities in the neighborhood of the 
shock are known, the configuration of the shock wave | 
itself is readily determined. In 
Returning to Fig. (2), the shock angle o is given by 2 | 
gam =~ UD, + a — Us — m) Fic. 4. Shape of ‘‘normal’’ shock wave in a flow involving a 
5 : . : roof-top velocity disturbance. 
or approximately 
. — ? , 
ae ; . dg Evaluation of Eq. (29) yields ”- 
o = (v — (/.6) (U; — U2) (27) . 
Pp 8M? A sin rmy cos mmy (30) 
The shock deviation from its original vertical directi — ; FFT . oo “seo gesag Se . 
8 lin cates Em? (Ui/U2) —1 _ m[l — (U2/U)] Th 
is seen to be 
Maximum values of P// occur at points given by 
co’ =a—5 (28) ae 3 - 
tan rmy = (U2/U;)(8/m)AM;? (31) for 
where Fig. 3 shows maximum values of P/F, or the ratio of - 
5 = hi" ) shock amplitude to wall amplitude, for various Mach 
6 = h'(—my). é ‘ oa er 
. Numbers, with y = 1.4. 1e dashed curve of thi 
‘ Numt tl 1.4. The dashed f this 
By integration, the equation for the shock deflec- S@me figure is a plot of maximum P/E with the effect lor 
tion P(v) becomes of vorticity neglected. It is interesting to notice that mP 
the effect of vorticity becomes large at high Mach 
. = 3 . € 
ss fay i Cros (29) Numbers, where the solid curve approaches a nonzero 
‘“ : r asymptote while the dashed curve tends toward zero. lor 
for 
1 
DISTURBANCE OF FINITE EXTENT 
Whereas the previous example was for the case of a disturbance extending to infinity, greatly simplified formulas 
can be derived if continuous disturbances of finite extent are considered. Since the entropy and vorticity are re- 
lated to 1 and its derivative, respectively, they can be written from Eqs. (20) and (21) as follows: lt 
' ” . pap 
(S — Ss)/R = Do(u;/ U4) (32) Nu 
Nut 
Q. = D,(Om,/Cn) 33) 
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where D,; and D2 are constants. Now if the characteristic waves intersect the axis x = 0 between the points Vo 
and yi, Eq. (13) becomes 
1 , U2 U2 OD: [ rs dé 
f= -tP — + Dz. — | u(y) + (y - nydra(n) f - is (34) 
2 U; U; 2 yo Je £2 + 9 9 


The integral term in Eq. (34) can be evaluated by integrating first with respect to and then by parts, provided 
that #(y) is continuous, as has already been assumed here, so that u:(vo) = m(¥1) = 0. The result is 


l CD, 32 U2 D U2 K 7 
g(v) = = — | my) = Amy) (30) 
2 8 2 "ee " 
say. 
Now Eq. (14) can be written, for x = 0, 

C.D “yt i édé K “1 uW4(n) 

2(0, Vv) = — — / duy(n) = = — ? = -dn (36) 
2x8 Jy £* +- (y — »)* 7B ~n Vu? 


Where the symbol ® indicates that Cauchy’s principal value is to be taken. The result is finite for continuous 
u(y), and after some simplification, one obtains: 


U2 A *¥ 4, (n) 
v2(0, vy) = 1 38M,? rf =. dn (37) 


1 7 


Example: Roof-Top Disturbance 


As an example, a roof-top velocity disturbance, first used by Howarth,’ will be considered. Let 


ui(n) = —(U;/m)(mn)e, 0 <7 < (1/m) 
u(y) = —(U,/m)(2 — mnje, (1/m) < 9 < (2/m) 
u(y) = 0, 0 > 7 > (2/m) 
In this case, 
“ U4(n) Ure 
ef : dy = — [my In my! — 2(my — 1) In |my — 1) + (my — 2) In my — 2] 
~~ =e m : 


The equation for the shock deflection P(y) is now determined by Eq. (29), and the constant of integration is ad 
justed so that Pis zeroat my = 1. The following constants are defined: 


G = ABM,?/rm[(U,/ U2) — 1], H = 1/[1 — (U2/U))] 
The equations for the shock are then found to be: 
mP/e = —G[—my In (—my) — 2(1 — my) In (1 — my) + (2 — my) In (2 — my) + 2] + (4/2) 
for my < 0; 


mP 


; (my)? l 
= —G[my In my — 2(1 — my) In (1 — my) + (2 — my) In (2 — my) + 2(1 — my)| — a| eee: 
€ 


for0 < my < 1; 


mP ; (my)* 3 
= —G[my In my — 2(my — 1) In (my — 1) + (2 — my) In (2 — my) + 2(my — 1)| — HM] 2my — i 


9 


forl < my < 2; and 
mP/e = —G|[my In my — 2(my — 1) In (my — 1) + (my — 2) In (my — 2) + 2] — (7/2) 


for my > 2. 
The quantity mP/eis plotted against my in Fig. 4 for three different Mach Numbers, with y = 1.4. 


CONCLUDING REMARKS For example, as already pointed out in connection with 

the sinusoidal case, the subsonic flow is rendered parallel 

It can be determined from the equations of this immediately behind the shock wave. This is true for 
paper that the particular free-stream supersonic Mach all continuous velocity profiles and would apparently 
Number V(y + 3)/2 has certain unique properties. occur also in a bounded stream. It results from the 
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opposing effects of the vorticity and source terms in the 
constant A.* 

At this same Mach Number, it also occurs that 
P(y) is proportional to h(x)—1.e., the shape assumed 
by the distorted shock wave is similar to that of the dis- 
turbance-producing wall. This case is included in 
Fig. 4, where the dashed line refers to M,= Vy + 3)/2 
= 1.48....(forair). The shape of the wave at 
this Mach Number is made up of two parabolic arcs, 
as is the wall that produces a roof-top disturbance 
profile. 

* The unique property of this Mach Number in the case of 
flow past a curved wall with a normal shock wave—namely, that 


the downstream radius of curvature is infinite—is apparent in the 
results of Emmons‘ and Tsien.® 
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The Instability of Pitching Oscillations of an 


Airfoil in Subsonic Incompressible 
Potential Flow 


BENJAMIN SMILG* 


Wright-Patterson Air Force Base 


SUMMARY 


A theoretical study is carried out to determine the conditions 
under which steady-state pitching oscillations of an airfoil in one 
degree of freedom are possible in two-dimensional subsonic in- 
compressible potential flow. The results show that such oscilla- 
tions can occur if the airfoil moment of inertia is sufficiently great 
and if the axis of rotation is located at a point that is forward of 
the airfoil quarter-chord but not too far forward of the airfoil 
The practical significance of these results with 


leading edge. 
and airplane stability problems is briefly ex- 


respect to flutter 
amined 


INTRODUCTION 


[' HAS COMMONLY BEEN ACCEPTED that, in order for 
flutter to occur, there must exist two or more de- 
grees of freedom that are coupled together by aerody- 
namic, elastic, or inertia effects. However, specialists 
in this subject have recognized for some time that the 
linearized aerodynamic theories for an oscillating air- 
foil in potential flow, both at subsonic incompressible 
and at supersonic speeds, indicate that undamped os- 
cillations of an airfoil in one degree of freedom (pitch) 
are possible under certain conditions.! In reference 2 
there is described the oscillation of ailerons in one de- 
gree of freedom at transonic air speeds, a phenomenon 
that has already been encountered in practice. Except 
for a brief study made in reference 3 to investigate 
single degree of freedom oscillations with respect to 
flutter, practically no work has been done on this sub- 
ject at subsonic incompressible speeds. It appears 
desirable to carry out further studies of this problem in 
order to define more clearly the conditions under which 
such oscillations may occur. In this manner it may be 
possible to explain why this phenomenon has not been 
encountered in practical aircraft and to show under 
what conditions it might occur. This is particularly 
advisable in view of the extremely unconventional 
aircraft configurations that are being considered at 
the present time in connection with high-speed 
flight. 

It is important to distinguish between oscillations of 
the type discussed in this report which occur in poten- 
tial airflow and those oscillations that occur due to 
separated airflow (e.g., wing stall). In both cases, the 
oscillations can be said to be caused by a time lag be- 
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tween the airfoil motion and the aerodynamic pressures 
acting on the airfoil. In the first case, this time lag 
can be predicted by potential airflow theory and occurs 
principally because of the effects on the aerodynamic 
circulation induced by the vorticity shed by the airfoil 
as it oscillates. However, in the latter case, the nature 
of the flow changes from relatively smooth stream lines 
to a separated type as the airfoil angle of attack is in- 
creased, the time lag being due primarily to the separa- 
tion phenomenon. 


DERIVATION OF EQUATIONS DEFINING CONDITIONS FOR 
STEADY-STATE OSCILLATIONS IN ONE DEGREE OF 
FREEDOM 


Assuming two-dimensional flow, the equation of 
equilibrium describing the steady-state pitching oscilla- 
tions of a rigid airfoil about a fixed axis of rotation is, 


according to reference 3, 


ig We F ° , 
ae o (1+ jg) |} + M,’ =0 (1) 
pb w 


where 

I, = mass moment of inertia of the airfoil about 
its pitching axis (slug-ft.? per ft. of span of 
the airfoil) 

p = air density (slugs per cu.ft.) 

6b = semichord of the airfoil in ft. 

®_ = circular natural frequency of the airfoil in 
pitch about the fixed axis of rotation at 
zero airspeed (rad. per sec.) 

w = circular frequency of the airfoil oscillations 
in pitch (rad. per sec.) 

g = structural damping coefficient 


j #=v-!1 

nondimensional complex coefficient for the 
aerodynamic pitching moment acting on an 
airfoil that is oscillating about a fixed axis. 
Positive values correspond to moments in 
the same direction as the displacement. 
The actual moment in ft. Ibs. per ft. of 
span is obtained by multiplying 7,’ by 
the term mpb‘w? (see reference 4) 


M,,’ can be described in terms of the coefficients de- 
fined and listed in reference 4 by the formula 
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For the purposes of this report, it will be assume 
that the control surface and tab deflections are zero~ 


FIG 1A 
REGION OF INSTABILITY IN ONE DEGREE OF FREEDOM 
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where 


a 


6 = 


location of the axis of rotation measured from 
the airfoil mid-chord point, as a fraction of 
the airfoil semichord (positive when located 


aft of the airfoil mid-chord point) 


angular displacement of the airfoil in pitch 


(positive when trailing edge moves down) 
angular displacement of 


aileron with respect 


to the airfoil (positive when trailing edge 


moves down) 
angular displacement of tab with respect 


to 


the aileron (positive when trailing edge 


moves down) 
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x}O 


FROM LEADING EDGE ( 








LOCATION OF AIRFOIL ROTATIONAL AXIS MEASURED 


500000 ~—~750,000 6 1p00ne0 


| I 
o 250,000 


AERODYNAMIC COEFFICIENT *-[M'x Ip 








assumed 
re zero— 
 simplif. 


le imagi- 
T, an ex. 
ws that, 
naginary 
y part of 
ty (v/bw 
r higher 
ore, the 


REEDOM 





HE 





PITCAMNG 


term (1/2) + a can be positive or negative depending 
on the rotational axis location. It is therefore to be 
expected that the imaginary part of 7,’ can become 
positive, corresponding to a condition of instability, for 
certain values of the reduced velocity and the rotational 
axis location. 

Eq. (1) can be rewritten in the following form for 
the case where the structural damping coefficient is 
zero: 


- w 
) 


fi -(2)] + tte +5 = 0 ©) 


where [1/,']g= real component of 1/7,’ (in phase with 
a) and [M,'], = imaginary component of VM,’ (90° 
out of phase with a) so that 


(I,,/mpb*) [1 — (wa/w)?] + [Ma’Je = 0 (4a) 
[M,']; = 0 (4b) 


It may be noted for the case #, = 0 that 
mpbt = —[M,'Jr (5) 


Graphs showing the relationship between the re- 
duced velocity v/bw and the location of the rotational 
axis required to satisfy Eq. (4b) are shown in Fig. 1 
for an airfoil with zero aileron and tab motions relative 
to the main airfoil. Points on the boundary corre- 
spond to an exact solution of Eq. (4b); points within 
the cross-hatched region correspond to a condition of 
dynamic instability. Fig. 2 is a plot of —[M,’Jr 
against the location of the rotational axis. The solid 


I 


a 


boundary line corresponds to the solid line of Fig. 1— 
i.e., the locus of points for which Eq. (4b) is uniquely 
The diagonal dotted lines correspond to 
In Fig. 3, the rela- 


satisfied. 
lines of constant reduced velocity. 
tionships shown in Eqs. (4a) and (4b) are plotted for 
constant values of the inertia parameter; the rota- 
tional axis location is used as the ordinate and a non- 
dimensional stiffness parameter w,b/v as the abscissa, 
where w, is the natural frequency (at zero air speed) of 
the airfoil about the rotational axis expressed in radians 
per second. Points that are located within the bound- 
aries of the curves correspond to conditions of dynamic 
instability. This figure is of the same form as that 
shown in reference 3, but minor discrepancies have 
been corrected and the inertia parameter has been ex- 
tended considerably to cover a range of values not 
only applicable to flutter but also to airplane and 
missile stability problems. The determination of 
whether a given system will be stable or unstable can 
easily be carried out by the use of Fig. 3. 


ESTIMATION OF OSCILLATION FREQUENCIES 


Eqs. (4a) and (4b), which define the conditions under 
which airfoil steady-state oscillations of constant ampli- 
tude will occur if no mechanical damping is present, are 


(I.,/mpb*) [1 — (we/w)?] + [Ma'lr = 0 
[M,'|; = 0 
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It is obvious that this set of equations can be solved 
for the values of the reduced velocity, inertia, and rota- 
tional axis location for which the oscillations will be 
marginally stable and the oscillation frequencies also 
determined. The results of such solutions are already 
shown in Figs. 1 through 3. However, the oscillation 
frequency can be more easily computed by an alter- 
native method based on replacing the real oscillatory 
aerodynamic moment by the static one. The aero- 
dynamic moment acting on an oscillating airfoil‘ is 


aerodynamic moment = mpb‘w*[M,’]a (6) 


If we consider only the case of constant amplitude os- 
cillations—i.e., [17,']; = 0, we can then write 


aerodynamic moment = mpb‘w?[M/,’]ra (7) 


It is easy to show for the static case (w = 0) that M7, = 
L, = Oand that Eq. (6) can be simplified into 


static aerodynamic moment = 27pv°b?[(1/2) + ala 
(8) 


The question then arises as to the error introduced by 
using the static aerodynamic moments of Eq. (8) in 
place of the theoretically correct real values of Eq. (7). 
This error can readily be seen by plotting the ratio of 


these moments. Thus, 
b*w? 
. ° [Ma'Ie 
real oscillatory moment v? (9) 
n = 7 . = - eed . 
static moment 2[(1/2) + a] 


This ratio, shown in Fig. 4, is close to unity for the range 
of parameters Consequently, the real 
moment for the oscillatory case can, in general, be ap- 
proximated by that for the static case over the range 
of pertinent reduced velocities (above 23.0). The 
maximum error is about 8 per cent at the lower values 
of the reduced velocity and becomes less than 2 per cent 
at values of the reduced velocity above 100. It is also 
interesting to note that the computations indicate the 
curve shown in Fig. 4 to be applicable regardless of 
whether the rotation axis is ahead of or behind the 
hinge line. 

The oscillation frequency w in radians per second can 
then be determined as follows: 


considered. 


elastic stiffness +- aerodynamic stiffness 


| 
= 4/ 
w \ i. 


| Dqy2 9 1 
«gy -"SS (10) 
b°(I,,/mpb*) 


It may be uoted that, if the rotation axis is ahead of the 
airfoil quarter-chord point, the term (1/2) + a@ will be 
negative, and therefore w will be greater than w,. The 
parameter n may generally be taken as unity, or it may 
be estimated from Fig. 4. A method of successive ap- 
proximations could be used, if desired, to obtain a more 


accurate value. A check of the oscillation frequencies 
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FIG. 3A 
REGION OF INSTABILITY IN ONE DEGREE OF FREEDOM 
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FIG. 3B 
REGION OF INSTABILITY IN ONE DEGREE OF FREEDOM 
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obtained from Eq. (10) with those given in Fig. 3 indi- 
cates excellent agreement. ° 
For the case w, = 0, Eq. (10) simplifies into 


wb/y = WV —2[(1/2) + a}n/(Iq/xpb*) (11) 


Eq. (10) thus defines the oscillation frequency pro- 
vided that the system is marginally stable dynamically. 
Therefore the values of a and J,/mpb‘ used in these two 
formulas must be chosen so as to correspond to a condi- 
tion of marginal stability—i.e., [17,’]; = 0 as shown in 
Figs. 1 through 3. 


ANALYSIS OF RESULTS 


The results of this theoretical analysis show that in 
order for undamped oscillations of the airfoil to exist 
at values of the reduced velocity greater than zero, the 
following conditions must be satisfied: 

(a) The axis of rotation of the airfoil must lie ahead 
of the airfoil quarter-chord but must not lie too far 
forward of the airfoil leading edge, the permissible 
amount depending on the magnitude of the reduced 
velocity (see Figs. 1 through 3). 

(b) The nondimensional inertia parameter J,/mpb* 
must exceed approximately 550. As the elastic stiffness 
restraining the airfoil from pitching or twisting is in- 
creased, the value of this inertia parameter must also be 
increased in order for these oscillations to occur. 

Study of these theoretical results show that the sub- 
ject single degree of freedom mode will not occur in 
most flutter problems where conventional airplane 
configurations are involved. Since the values of the 
inertia patameter J,,/mpb‘ for a full-scale wing, stabilizer, 
or fin are not likely to exceed about 30 at sea level 
and, besides, since the elastic axes usually lie aft of the 
airfoil quarter-chord, it would appear impossible for 
these airfoils to flutter in the torsional single degree of 
freedom regardless of the torsional rigidity of the 
structure. Even for a freely floating stabilizer or con- 
trol surface with the rotational axis located ahead of 
the quarter-chord, the inertia parameter of the full- 
scale airfoil would be too low to permit this mode of 
flutter, except perhaps in special cases where the rota- 
tional axis was at least a chord length ahead of the air- 
foil leading edge and where the airfoil was extremely 
heavy. However, the minijnum conditions required 
for flutter of an airplane structural component in this 
mode might possibly be satisfied in the case of a swept- 
back wing with an extremely heavy weight near the 
tip, an airplane empennage where the rear fuselage 
structure is flexible, and the ‘‘floating”’ airfoil of a small 
but relatively heavy wind-tunnel model or missile. 
In an actual flutter analysis for these cases, the effects 
of the other degrees of freedom including the rigid body 
modes would require consideration. 

The theoretical results of this report can also be ap- 
plied to problems involving airplane stability by con- 
sidering the elastic restraint at zero air speed to be equal 
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to zero. If the complete airplane is considered as a 
rigid body free to oscillate in one degree of freedom 
about an axis through or close to its center of gravity, 
the values of the nondimensional inertia parameter 
I,,/mpb* can easily exceed 550. A check of a number of 
typical airplanes show the following range of values for 
I,,/mpb4 at sea level, where /, is the mass moment of 
inertia of the entire airplane about the airplane c.g. 
per unit span of the airfoil in question: (1) for the 
vertical tail only, with respect to yaw-——2,000 to 
20,000; (b) for the horizontal tail only, with respect 
to pitch—3,500 to 18,000; (c) for the wing only, with 
respect to pitch—200 to 500. These values will 
increase at the higher altitudes, becoming approximately 
ten times as large as those listed above at an altitude 
of 60,000 ft. However, in order for the instability to 
occur, there must be an effective axis of rotation of the 
airplane which lies forward of the airfoil quarter- 
chord but not too far forward of the airfoil leading 
edge, the permissible amount increasing as the value of 
the nondimensional inertia parameter is increased. 
Even if this parameter were equal to 1,000,000, which 
is far above the value expected in any practical case, 
one degree of freedom oscillations would not occur if 
the axis of rotation was located more than two airfoil 
chord lengths ahead of the airfoil leading edge. With 
normal configurations, the airplane c.g. will lie far 
forward of this limit insofar as the tail surfaces are 
concerned so that this type of airplane instability due 
to the tail is not likely to occur. It is also to be 
expected with normal configurations that the airplane 
c.g. will be aft of the wing quarter-chord; thus, insofar 
as the wing itself is concerned, this mode of instability 
is not likely to be present. However, it can be visual- 
ized that some difficulty might be encountered with 
special configurations, involving sweptback wings or 
an extremely short tail length, and also with flexibly 
supported wind-tunnel models. It must be remem- 
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bered that this analysis is based on the assumption 
that the motion of the airplane can be described as 
pitching or yawing about an effective rotational axis 
passing through the airplane center of gravity. Before 
drawing any broad conclusions, it would be desirable 
to check the validity of the assumption by studying 
the oscillation modes of typical aircraft considering the 
coupling of pitch with vertical translation, the coupling 
translation, and roll, and the 


of yaw, lateral 
between the body 


more general coupling 
modes and the modes involving structural deforma- 


rigid 


tions. 

It is interesting to observe from Fig. 3 that an airfoil 
that is stable in this one degree of freedom at a given 
air speed may become dynamically unstable in the 
same mode at a higher air speed. This result, however, 
does not apply when the elastic restraint is zero 
(w, = 0). For this case, an airfoil that is stable in this 
mode at any air speed will theoretically be stable at all 
air speeds; on the other hand, if it is unstable at any 
air speed, it will theoretically be unstable at all air 
speeds and the oscillation frequency will be directly 
proportional to the air speed. 

The results given in this report are directly applicable 
only to the simple airfoil of infinite aspect ratio in 
subsonic incompressible flow and for a system with no 
damping except that contributed by the aerodynamic 
pressures acting on the airfoil itself. The extensions 
of the methods described in this report to cover other 
conditions are reasonably straightforward, provided 
that suitable expressions for the aerodynamic pressures 
can be determined. It is important to note that the 
pseudo-static* aerodynamic forces and moments that 
are usually used in airplane stability computations 
would not indicate any instability of the type discussed 
in this report. This would indicate that, in the general 
case, pseudo-static aerodynamic terms should not be 
used in airplane dynamic stability analyses if it is 
desired to obtain accurate results regarding the rate 
at which oscillations are damped. 

*“Pseudo-static” is taken to mean the static aerodynamic 
values corresponding to the instantaneous angle of attack, thus 
neglecting effects of the lag in aerodynamic circulation 
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CONCLUSIONS 


(1) Undamped oscillations of an airfoil in pitch 
involving one degree of freedom may occur in subsonic 
imcompressible potential flow, provided that the magni- 
tude of the mass moment of inertia of the airfoil and the 
location of the axis of rotation satisfy the conditions 
shown in Fig. 3 of this report. The existence of this 
type of instability would not be predicted by the 
pseudo-static aerodynamic forces now commonly em- 
ployed in airplane stability analyses. 

(2) The application of these results to practical 
cases requires an extension of this study in order to 
examine the effects of aspect ratio, compressibility, and 
the couplings with the other degrees of freedom, in- 
cluding modes involving rigid body motions as well as 
those involving structural deformations. However, 
based on the simplified study made herein, it appears 
that: 

(a) Oscillations of this type are not likely to occur 
in flutter problems except in special cases involving: 
(1) a sweptback airfoil with an extremely heavy weight 
near the tip, (2) an airplane empennage where the rear 
fuselage structure is extremely flexible, and (3) a 
‘floating’ airfoil of a small but relatively heavy wind 
tunnel model or missile. 

(b) Undamped oscillations of this type are not 
likely to be important in stability problems for air 
planes of normal configurations but may be encoun 
tered in: (1) sweptback wings, (2) designs where an 
extremely short tail length is used, and (3) flexibly 
supported wind-tunnel models. 
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Mixing of Any Number of Streams in a 
Duct of Constant Cross-Sectional Area’ 


ROGER WEATHERSTON?t 
Cornell Aeronautical Laboratory, Inc. 


SUMMARY 


A one-dimensional analysis is applied to the mixing of any 
number of streams of a compressible gas in a duct of constant 
cross-sectional area. By the proper selection of significant pa- 
rameters and by the use of convenient functions, the solution of 


the problem becomes simple and brief. 


ASSUMPTIONS AND NOTATION 


Assumptions 

1) The mixing takes place in a duct of constant cross-sectional 

area and is completed within this duct. 

2) Wall friction is neglected 

3) The fluid obeys the equation of state of a perfect gas. 

4) y and R are assumed to be the same for all merging flows, 
Symbols 

M = Mach Number 

D = M/{1+ [(y -— 

G = (1+ yM2)/{1 + [(y — 1/2] Me2}V/0-D 


2).2} +1)/2(7 —1) 


D , ; 
N= GT MV ULF [ly — D244 / + yA) 
7 
p = static pressure 
P = total pressure 
T = total temperature 
V = velocity 
A = cross-sectional area 
m = mass flow 
y = ratio of specific heats 
R = gas constant 
Subscripts 
1,2,3,..., ” = condition of streams 1, 2, 3, , natinception 


of mixing 
f = station at which mixing is complete 


INTRODUCTION 


I ONE-DIMENSIONAL FLOW, the condition of the 
stream at any section of a duct is completely de- 
fined by four independent parameters (e.g., two state 
parameters, the cross-sectional area, and the velocity). 
In dealing with the problem of mixing, the author has 
found it convenient to describe both the merging and 
final flows in terms of total temperature, mass flow, 
Mach Number, and cross-sectional area. At the sta- 
tion where mixing is complete, the first two of these 
parameters are simple additive functions of the corre- 
sponding quantities for the individual streams at 
the inception of mixing; the third parameter is ob- 

Received March 17, 1949. 
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. 


tained through the solution of a simple equation for the 
end Mach Number in terms of the known initial total 
temperatures, mass flows, and Mach Numbers of the 
merging flows; and the fourth one either is given, or 
may be obtained, as the sum of the stream areas of the 
merging flows. It has been found (e.g., references 2, 3) 
that the selection of any other set of state and flow 
parameters to describe the initial and final flows would 
lead to equations that could not easily be solved, even 
in the simple case of the “‘ejector’’ problem (mixing of 
only two flows), without the use of simplifying assump- 
tions that would restrict their validity. 


ANALYSIS OF THE MIXING OF 1 STREAMS IN A CON- 
STANT AREA Duct 


The analytical treatment of this problem is consider- 
ably simplified by the use of the Mach Number func- 
tions D, G, and N, which were introduced in reference 


1. A plot of these functions is shown in Fig. 1. With 
these functions 
m = Vy/R(PAD/VT) (1) 
and 
pA + mV = PAG (2) 


The momentum equation may be written in the form 
i=n 
) PAG; as PA Gy (3) 
i=1 


From Eqs. (1) and (3) one obtains 
i=n . 
> (mv T, ‘N,) = mW T,/ Ny (4) 


From Eq. (4), together with the continuity equation 


i=n 


m, = >> m 
i=1 
and the energy equation 
i=n 
mT, = >, m7 (5) 
i=1 


one obtains 
=n =n ‘=n J 
Ny = &>> m, > mT; /> (mV T,/N:i) (6) 
i=l §=1 / i=l 


Then M, can be obtained from Fig. 1 or from a tabula- 
tion of the function JN. 


697 








698 JOURNAL OF 


THE AERONAUTICAL SCIENCES-NOVEMBER, 


1949 





13 


09 


0.8 


0.7 


0.6 + 


0,G, N 


0.5 


0.4 F 


0.3}+-—+ 


0.2 


0.1 




















0 0.2 0.4 0.6 0.8 LO 12 


16 18 20 2.2 2.4 2.6 2.8 3.0 


MACH NUMBER 


EXAMPLE 


For the sake of illustration, this procedure will now 
be applied to the solution of an “‘ejector’’ problem. 


Both inducing and induced gases are air with y = 1.4 
and R = 1,715 ft.lbs. per slug R°. 
Given 
m, = 1 slug persec; 7; = 500°R.; JA, = 0.2 
m, = 0.5 slug per sec.; 72 = 1,500°R.; M: = 


find M,. 


Mach Number functions D, G, and N. 


Solution 

From Fig. 1, with the given values of AJ, and Mz, one 
obtains 

N, = 0.190, Ne = 0.456 
Eq. (6) gives 
V (1.5) [(1) (500) + (0.5) (1,500)] 
V 500/0.190 + 0.5V 1,500 0.456 
Therefore (Fig. 1), WM, = 0.303. 
(Continued on page 704) 


N, = = 0.271 
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RIEF REPORTS of investigations in the aeronautical 
B sciences and discussions of papers published in the JOURNAL 
will be presented in this special department. The publication will 
be completed 6 to 8 weeks after receipt of the material. No proof 
will be sent to the authors. The Editorial Committee does not hold 
itself responsible for the opinions expressed by the correspondents. 

Contributions should not exceed 800 words in length. 


Structural Damping 


W. Pinsker 
145, Cove Road, Farnborough, Hants., England 
July 27, 1949 


Ww REFERENCE TO the article by Walter W. Soroka with 
the title ‘‘Note on Relations Between Viscous and Struc- 
tural Damping Coefficients’? (JoURNAL OF THE AERONAUTICAL 
ScrENCES, Vol. 16, No. 7, p. 409, July, 1949), I should like to make 
the following comments. 

The physical meaning of the complex spring stiffness k(1 + ig) 
as used by Mr. Soroka in his paper on structural damping has been 
misinterpreted: a correct interpretation may assist in the explana- 
tion of the contradictory conclusions of the paper. The calcu- 
lated increase in frequency arising from the introduction of the 
term &(1 + zg) does not result from the reciprocal dependence 
of the damping coefficient on the frequency w; such a relationship 
would only cause a further decrease in the natural frequency of 
such a system. The imaginary part of the complex stiffness 
ikg does not generally represent a force in phase with the velocity 
which in the case of a damped oscillation is 

[(7/2) + tan-! (6/27)} 
ahead of displacement. 

By splitting up &(1 + ig)x into the actual components in phase 
with displacement and velocity, ‘‘effective’’ restoring- and damp- 
ing-coeflicients can be obtained; these are, respectively: 


»\2 2 > 
fie) 
w) 2m Ox 


c’ = (k/w)g = OP/dx 


Substituting these quantities into the normal differential equation 
for a system with viscous damping 


méi+c’zr + k’x =0 


gives the same results as obtained in the original paper: 


ee | 


Ihe ’ - cauieesaanediiniieaia 
wail |k ee k / — 
V m 4m? V 2m -? v1 wi 


b= Ong + Vite 


Now the increase of the natural frequency with g can simply be 
explained by an increase of the effective k’ with g; the increase 
in natural frequency due to this turns out to be just twice as much 
as the normal reduction in frequency due to the damping term c’. 

It appears that the complex structural stiffness used has been 
evaluated from steady oscillations. Further information on the 
actual magnitude and phase of the complex spring stiffness during 
a damped oscillation of a system with structutal damping has to 
be known before drawing any conclusions on the nature of damp- 
ing term to be used in the differential equation. 


On the other hand, observation of the natural frequency of a 
system with variable structural damping will also suffice to deter- 
mine the actual working of the damping component. If, for 
instance, a damping coefficient c’ = kg/w would represent com- 
pletely the structural damping effect without changing the spring 
stiffness k, the frequency must be expected to give 


- Jf ” = 
wn = Vk/2m V 1471 — 32 


Heat Addition in a Pipe 


D. G. Samaras 
July 7, 1949 


a A PAPER “On the Addition of Heat to a Gas Flowing in a Pipe 
at Subsonic Speed’”’ (JOURNAL OF THE AERONAUTICAL ScI- 
ENCES, Vol. 16, February, 1949), Drs. Foa and Rudinger treated a 
specific case of a generally well-known and understood problem on 
the basis of the following erroneous assumption: 

**.,. Contrary to widely accepted opinions: 

‘“‘(a) The addition of heat to a subsonic flow modifies the flow 
condition both downstream and upstream of the region of heat 


” 


addition. ..,’’ ete. 

To the writer’s knowledge, such misconceptions are absent at 
least in his own paper! (quoted by the authors). On the con- 
trary, the misunderstanding seems to be due rather to a misinter- 
pretation of the literature by the authors. 

The present writer has treated in integrated form the general 
problem, not the simplest case considered by the authors, taking 
into account frictional losses, general configurations of converg- 
ing and diverging ducts, and subsonic, as well as supersonic, 
flows. 

This problem arose in 1943 when the writer was entrusted with 
the investigation of the best power plant for the M-52 supersonic 
aircraft. Up to that time the matching of turbojets and ducted 
fans was performed by disregarding the characteristics of the 
combustor and slightly modifying the compressor characteristics 
to take into consideration the pressure losses in the combustor. 
Such a process was insufficient for matching the configurations 
considered at that time® * —namely, turbojets and ducted fans 
with afterburning and ram-jets. 

The problem consisted in finding the variation of the total head 
and static temperature and pressure ratio as a function of the in- 
let and exit conditions, the rate of heat addition, and the geom- 
etry of the duct. The latter two were combined in one area 
parameter m, and the following functions were obtained: 

To 
= F,(M,, Mo, m), 


lt 2t 


F,( M,, M2, m) 


Analogous relations were found for the static ratios. 

These functions represent a family of surfaces in a three-dimen- 
sional space; to plot them in a two-dimensional graph, current 
engineering practice suggests a set of constant values of the one 
variable—in our case m—and in each of these graphs a series of 
values for one of the remaining variables—in our case, the inlet 
Another choice could be to keep constant the 
This isa 





Mach Number. 
exit Mach Number or temperature and pressure ratio. 
matter of convenience and taste. 

Fig. 1 shows the combustor characteristics evolved at that time, 
which is nothing but a plot using slightly different coordinates, 
although suitable for the matching procedure, from those of 
reference 1. From this and Figs. 3, 5, and 7 of reference 1, it is 
obvious that the combustor inlet and exit boundary conditions 
must vary to fit the heat addition process and, in addition, the 
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SIMULATION OF THE INTERNAL AND EXTERNAL FLOW /N 
A COLD PRAMIET BY MASS ADDITION. 











flame stability limits (not even mentioned in the authors’ paper) 
The same figures show the distribution of the static and total head 
temperature and pressure along the duct at certain inlet and exit 
conditions corresponding to a given heat addition and may also be 
used to find the inlet conditions when the exit conditions and 
temperature rise are given 

Later, in the summer of 1944, when the fragments of the ex 
perimental V-2 were brought from Sweden, use of these results 
was made in order to investigate the heat-exchange processes 
taking place in the rocket nozzle (supersonic heat addition and 
extraction—a case that has not even been touched by the 
authors!). 

As in the above project ram-jets were investigated,” the writer 
suggested an experimental arrangement to simulate within a first- 
order approximation the internal and external flow in a ram-jet to 
be tested in a cold tunnel (Fig. 2). The general principle used 
was that of substitution of mass addition for heat addition. 

It should be mentioned here that early during the war experi- 
mental evidence showed that when choking occurs in the turbine 
nozzles of a turbojet, the variation of the temperature rise in the 
combustor is always accompanied by a change of the inlet condi- 
tions of the combustor. This process is obvious from my paper,! 
and the physical explanation was given lately by Kahane and 


Lees. 
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Authors’ Reply 

Dr. Samaras’ opinion (above), that modifications of upstream 
conditions in a heated flow are peculiarly associated with the 
occurrence of ‘“‘choking,’”’ is indeed at variance with our conclu- 
sion! that such modifications will a/ways occur when the heating 
rate is changed in a subsonic pipe flow. His introductory state 
ment can hardly be reconciled with the undeniable existence of 
this divergence of views. 

Our paper dealt with the subsonic case exclusively. The super 
sonic case, which is entirely different, was treated in a separate 
publication. 

A full answer to the remainder of Dr. Samaras’ comments will 
be found in our recently published reply to other comments 
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Comparison of the Rayleigh-Ritz and 
Complementary Energy Methods in 
Vibration Analysis 


Paul A. Libby and Robert C. Sauer 

Department of Aeronautical Engineering and Applied Mechanics, 
Polytechnic Institute of Brooklyn, Brooklyn, N.Y. 

August 29, 1949 


E REISSNER’S CONTRIBUTION to the Readers’ Forum! sug 
¢ gested that the results of an investigation of the relative 
accuracy and computational difficulties of the Rayleigh-Ritz and 
complementary energy methods applied to the lateral vibration 
of a beam might be of some interest. The comparison is made for 
the symmetric modes of vibration of a free-free beam with a 
concentrated mass at its centerline of symmetry. In this note 
the stationary conditions to be satisfied approximately by the 
assumed modes of deflection are formulated from the differential 
equations of motion for the general case of a variable section beam 
The results of the two methods with several assumed mode shapes 
for the fundamental and second modes of vibration are then com 
pared for the case of a uniform beam with the exact solution 

With the assumption of harmonic motion of circular frequency 
w, the ordinary differential equation describing the modes of 
lateral deflection Y of a beam is 


(EIX")" — pwtX = 0 (1) 


where E/ is the bending rigidity, p is the mass density per unit 
length of span, and the primes denote differentiation with respect 
to the space variable along the beam, x. If only symmetric modes 
of vibration are considered and if the origin of the x axis is placed 
at the centerline of symmetry, then only one-half of the beam 
(0 < x < L, where L is the semispan) needs to be considered 





Th 
bu 


(Fe 


Eq 
apy 
fror 
tha 


The 
app 
vib1 
Stat 

if 


tion 


whe 
The 


w[J~ 


and 


Jour. of 
u person 


Ram-Jet 
*, 1946. 
t Aircrafi 


stream 
ith the 
conclu- 
heating 
y state 


ence of 


> super 


eparate 


nts will 
its 


Flowing 
Vol. 16, 


Flowing 
port No 
9 

Readers’ 


6-567 


hanics, 


sug 
lative 
tz and 
ration 
de for 
vith a 
+ note 
yy the 
ential 
beam 
hapes 
} com 
yn 

iency 
les of 


(1) 


- unit 
spect 
nodes 
laced 
beam 
ered 


AS.M.E 


approximately by the usual Rayleigh-Ritz method. 
from Eq. (1 
that 


Thus, Eq. (3a) may be 


7 < 


READERS’ FORUM 701 
TABLE 1 
Exact Rayleigh-Ritz Method Complementary-Energy Method 
Mode shape a | Mode shape b | Mode shape c | Mode shape a | Mode shape 6 Mode shape « 
} 
M/pL a an a TA a 7 & an %A a % & a "oA a cA a A a A 
0 5.60 | 30.25 | 5.72 2.14 5.60 0 31.28 3.41 | 6.70 19.2 5.67 1.25 5.60 0 30.85 1.98 | 5.61 0.17 
4.56 24.50 4.61 1.09 4.56 0 27 .86 13.70 >. 61 23.0 4.60 0.88 4.56 0 25.00 2.04 | 4.57 0.22 
1 4.12 22.84 4.25 3.27 4.21 2.31 26.20 14.50 5.25 27.4 4.23 2.79 4.20 1.82 23.70 3.76 4.24 3.02 
) 3.80 | 22.65 3.96 :.2 3.94 3.69 25.16 11.0 4.95 30.2 3.95 3.95 3.94 3.69 23.10 1.98 3.95 3.95 
3.52 22.40 3.53 0.50 3.52 0 23.90 6.70 | 4.46 26.7 3.52 0.17 3.52 0.28 22.70 1.34 3.53 0.50 
Relative work 1.00 41.00 4.00 0.25 1.50 6.00 6.00 0.40 
, 9 . lary p aw 2 “J pee | 
Ifa concentrated mass 2./ is located at x = 0), then the boundary Se ; 
24s > aj y\ w* (1/E/) dn pe( edt x 
conditions are j ) 
j 1 ; 0 oJ x n 
XO) = XL) = XL) = O (2a) le . . 
dn per(E)dé | dx pe en dx + 
(EILX"’)’x-0 = w?MX(0) (2b) J % ” 7! 
{| 
The last boundary condition may be expressed in an equivalent Med Oreu(0) | ( = i,k L2z p (5b) 


but more convenient form as 
°L 

wVX(0) = —w? pX dx (2c) 

(For a reference on these boundary conditions cf. M. Goland and 
Y. L. Luke, ‘“‘The Flutter of a Uniform Wing with Tip Weights,” 
Paper No. 47-A 
Multiplying Eq. (1) by 6X, 


36.) 


integrating by parts, and consider 


ing Eqs. (2a)-(2c) with respect to both Y and 6X, one obtains 
the stationary condition 
| «> 
54 (1/2) EI X’'*dx (w?/2) X 
| 0 
“E 
pX*dx + MX*(0) = 0 (8a) 
0 
Eq. (3a) expresses the stationary condition that is satisfied 


However, 


) with the boundary conditions [Eq. (2a)] one finds 


yw S "| 
El) | dn | pX(é&)dt 
A n 


written with Eq. (4) as 


x” = 


a | a» 3 L 2 
(1/EI) | dn f pX(t)dt | dx — 
0 x ” 


(1/2) pX*dx + MX*(0) ' = 0 (5 


0 


Comparison of Eqs. (8b) and (5b) indicates the increased ac 
curacy that, as Reissner pointed out, could be expected from the 
complementary energy method; the assumed mode shapes are 
integrated three times, while in the Rayleigh-Ritz method they 
are differentiated twice and integrated once. On the other hand, 
if, as is frequently done for simplicity, polynomials are used for 
the ¢ functions, the computational work is greater with the com 
plementary energy method, since the terms in the ¢’s do not drop 
out upon integration as they do with differentiation. It was for 
the purpose of comparing the relative accuracy and computa 
tional work involved in the two methods that this investigation 
was undertaken 

For the purpose of comparison the numerical work was carried 
out for a uniform beam for which the exact frequencies may be 
found by direct integration from Eqs. (1) and (2 In both the 
the natural frequencies were 
found in terms of the nondimensional ratio J /pl. 
» © the frequencies 


exact and approximate solutions, 
This provided 
a check on the exact solution, since as //pl 
reduced to those of a cantilever beam of total mass p L and for 
W/L free-free beam of 


semispan L and of total mass 2pL. 


= ( to those of a completely uniform, 


In the approximate solutions the mode shapes were assumed 
as simple polynomials of the following forms: 
Mode shape a: 


stationary condition presented by 





The stationary condition given by Eq. (5a) forms the basis of the 
application of the complementary energy method to the lateral 
vibration of beams and corresponds for this application to the 


E. Reissner. 
In both methods approximate solutions to the stationary cordi 


tions given above are obtained by assuming that 


p 


pe aje;(x) 


j=1 


X = 


where each ¢,(x) satisfies the boundary conditions on X [Eq. (2)]. 
Then Eqs. (3a) and ( 


5a) become 


Ee 


I 
Af E T¢;"" ox’ ‘dx — w? PLIF: dx + 
= 0 


Me,(0)¢x(0) (3b) 


—_—— 


and 


XY = a,[t' — 4&5 + 6: 6/A5(1 + 8) 
Mode shape pb: 
X = alt — 483 + 62 — 6/5(1 + B)] + aol — 1083 + 
uF: 13/3(1 + 8B) 
Mode shape c: 
XY = a(t? — 1/3(1 + 8) 
where — = x/Land B = M/pL. 


Mode shapes a and 8 satisfy all of the boundary conditions of 


equations (2a — c); mode shape c satisfies only the boundary 


conditions X’(0)=0,.M/ X(C) + pXdx = Oand since ¥’’’(x) 
0 
Howevet, it was considered to be of interest 


san be obtained with the two 


=0,X''(L) =0 
to determine the accuracy that 
methods with simple assumed mode shapes not satisfying all of 
the boundary conditions. The results are shown in Table 1. 
They are expressed in terms of the pure numeric a; where 


(a;/L?) WEI /p 


and where the 7 denotes the jth mode 

Examination of Table 1 indicates the increased accuracy of the 
complementary energy method. Especially noteworthy is the 
accuracy obtained by this method with mode shape c and in the 
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second mode with mode shape b. It will be noted that with only 
two ¢ functions—i.e., the minimum number—sufficient accuracy 
for engineering purposes is obtained for the second mode. To 
obtain similar accuracy with Rayleigh-Ritz at least three ¢ 
functions must be used. 

Also listed in Table 1 is an estimate of the relative work required 
for each calculation as compared to that for the conventional 
Rayleigh-Ritz calculation of the first mode (mode shape a). 
It will be noted that the work with the complementary energy 
method is approximately 59 per cent greater than that for the 
Rayleigh-Ritz method. However, for simplified mode shapes 
(mode c) and for the determination of the higher modes (mode 3, 
for example), this additional work is warranted by the increased 
accuracy. 

The results given here seem to indicate the superiority of the 
complementary energy method over the Rayleigh-Ritz method 
for certain vibration problems and are in agreement with the 
suggestions of E. Reissner regarding the application of the com- 
plementary energy method to flutter calculations. 
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A Simple Method of Calculating the Subsonic 
Rolling Characteristics of Sweptback Wings 


Edward N. Spryer 

Unit Head, Loads Group, Airplane Division, Curtiss- Wright Corpora- 
tion, Columbus, Ohio 

August 11, 1949 


ge FOLLOWING FORMULAS have been derived! for the rapid 

determination of the damping in roll parameter C;, and the 

aileron rolling effectiveness Ci; of sweptback wings at subsonic 
a 





speeds: 
C oc; 114.6 Aap cos Azg 
 ieadiael = - >? 22 AC 
: pb Sb AE,’ + 36.48a0 cos Ayz 
age 
2) 
~ 2y 
cr— ay 2) 
0 b 
. oc; Aaecos Arg 2 
Cig ee ee * (a; cos Ay,) = 
a 05,4 A Ez. aa 36.4840 cos Arg “a Sb 
Yo 
cydy (2) 
yi 
(Ci, )ae AE,’ + 36.48ae cos Azz 
Oe ——— (3) 
(Cp) = 0 AE,'V/1 — M* cos? Azg + 36.4809 cos Azg 
where 
C; = rolling mroment coefficient, L/gSb 
L = rolling moment, ft.lbs. 
q = dynamic pressure, lbs. per sq.ft. 
S = wing area, sqft. 
b = wing span measured perpendicular to plane of symmetry 
ft. 
p = rolling velocity, rad. per sec. 
\ = true air speed, ft. per sec. 
{ = wing aspect ratio, (b?/S) 
a9 = section lift curve slope, per deg. 
I,’ = effective edge—velocity correction factor for rolling 


moment 
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( = wing chord measured parallel to plane of symmetry, ft 

y = perpendicular distance from plane of symmetry, ft 

Ml = free-stream Mach Number 

A = angle of sweepback, deg. 

5, = aileron deflection measured in plane perpendicular to 
aileron hinge axis, deg. 

a; = rate of change of section angle of attack with aileron de- 

a 


flection, (Oa 06,) 
and the subscripts are: 


LE = wing leading edge 

HL = aileron hinge line 

t = inboard end of aileron 
0 = outboard end of aileron 


The above formulas are based upon the following assumptions 
(1) The airfoil section is constant along the span; (2) the right 
and left aileron deflections are of equal and opposite magnitude: 
and (3) the Prandtl-Glauert rule applies up to the critical Mach 
Number. 

In all three formulas the term AE,’ appears. 
this value with aspect ratio A appears in Fig. 22 of reference 2, 
It should be noted that in Eq. (2) the aileron deflection 6, is only 
that of one aileron; however, the resulting rolling moment coeffi- 
The parameter 


A variation of 


cient accounts for the effect of both ailerons. 


a, is a two-dimensional value and is dependent upon aileron- 
a 


chord ratio, trailing-edge angle, and balance and seal characteris- 

tics. No formula is given for the effect of Mach Number on 

Ci, A study of existing wind-tunnel data for sweptback 
a 


wings indicates only a slight variation of C1; with Mach Number 
a 


at Mach Numbers below the critical. A probable explanation 
for this phenomenon is that ap increases with increasing Mach 
Number while a, decreases. 
a 


A comparison with existing wind-tunnel data indicates that the 
theoretical values of Ci, and Ci, are slightly higher (within 10 
a 


per cent) than the experimental values. However, it is felt that 
in the preliminary design stages of an airplane (i.e., before wind- 
tunnel results are available) the above formulas will predict 
fairly accurate subsonic rolling characteristics when (1) the wing 
aspect ratio is greater than 3, (2) the sweepback of the leading 
edge is less than 50°, and (3) the Mach Number is below the 


critical. 
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Errata: On the Oscillating Rectangular 
Airfoil at Supersonic Speeds 


John W. Miles 

Department of Engineering, University of California at Los 
Angeles 

July 22, 1949 


_— FOLLOWING ERRORS appeared in my Readers’ Forum 
item in the June issue of the JOURNAL OF THE AERONAUTICAL 
SCIENCES (Vol. 16, No. 6, p. 381): 

The exponent in Eq. (2) should be —ixMx and not —2ix 


ra) 
_ ( tke t) (6) 
Ox 


if 


x 


Eq. (6) should have read 


a(x, t) = 
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FIG. 1, TYPICAL VARIATION WITH MACH NUMBER OF 
VELOCITY AT ANY POINT ON A STREAMLINE 
BODY OF REVOLUTION AT ZERO INCIDENCE 
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On the Compressibility Correction Factor for 
Axially Symmetric Bodies 


A. D. Young 
Senior Lecturer, Department of Aerodynamics, College of 
Aeronautics, Cranfield, England 


August 12, 1949 


+ poe CONTRIBUTION under the above title by W. H. Dorrance, 
appearing in the Readers’ Forum of the July, 1949, issue of 
the JOURNAL contains, in my view, two fundamental errors. 
Firstly, the potential function and its derivatives for an axisym- 
metric body becomes infinite on the axis within the body, and one 
cannot replace the values of the derivatives on the body by the 
values on the axis, as one can in two-dimensional flow. Hence 
Eqs. (7) and (14) in Mr. Dorrance’s note need to be modified to 
have any meaning; the integral in the latter equation as it stands 
is infinite. Secondly, toward the end of the note there appears 
to be some confusion between the incompressible flow on the 
body and that on the slimmer body defined by the transformation 
of Eq. (3) of his note. 

Eq. (12) is correct as far as it goes, and it tells us that Af(x) 
is independent of Mach Number for flow about a given body. 
One can deduce from this that the equivalent source-sink dis- 
tribution of a body is independent of Mach Number, a deduction 
that has already been made by Lees! and is also given in a paper 
of which the present writer was part-author.? To determine the 
pressure distribution in compressible flow about a given body, 
however, we must first determine the relation between that pres- 
sure distribution and the distribution at zero Mach Number on 
the slimmer body defined by the transformation of Eq. (3). 
Adopting Mr. Dorrance’s notation (i.e., using dashes to denote 
quantities associated with the latter body) we have tan @ = (1/8) 
tan 6’, and, hence, 


(b¢/U)y +» « = (1/8)(6'r"/ Uy? > 


But, from Eqs. (12) and (3) of Mr. Dorrance’s paper 


FORUM 703 


(dr)y +e = Kf(x:)/2xr = Kf(x1)8/2xr’ 
('r")y? —>,. = f(x1)/2xr’ 
It follows that K = 1/8?. 


Thus, the perturbation velocity in the axial direction on a thin 
body in compressible flow is 1/8? times that in incompressible 
flow on the thinner body obtained by reducing the lateral dimen- 
sions of the original body in the ratio 8:1. I understand that this 
result was first obtained by Géthert. 

This rule has been applied to determine quantitatively the 
effect of compressibility on the pressure distributions on typical 
streamlined bodies of revolution, and the results are sufficiently 
consistent to be summarized as shown in Fig. 1 herewith.? It 
will be seen that, though the effect is somewhat less than that 
given by the Glauert law, it is by no means negligible in com- 
parison. 
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Errata 


Theodore R. Goodman 

Associate Research Aerodynamicist, Cornell Aeronautical, Inc., 
Buffalo, N.Y. 

August 12, 1949 


ls MY PAPER “The Lift Distribution on Conical and Nonconical 
Flow Regions of Thin Finite Wings in a Supersonic Stream”’ 
(JOURNAL OF THE AERONAUTICAL SCIENCES, Vol. 16, No. 6, June, 
1949), I have noticed that the derivation of Eq. (10), and other 
equations [Eqs. (34) and (43)] of which Eq. (10) is the prototype, 
is incorrect. Eq. (10) itself and all subsequent results are cor- 
rect. I shall herewith derive this equation correctly using the 
notation given in my paper. 
Eq. (2) reads 


up? Z AC,(X1, V1) dX, dV, _ @) 
U 4n | , (X —X)? -— XY — ¥,)? — z73/72 


Differentiating with respect to Z, 





au/U) et | (ff — aG(X, Yi) dXi dX, | 2 
oZ 4 , (X — Xi)? — BY — Y,)? — 6927]3/ 





, B ‘ Z? ACP(N, Vi) dN d¥, 4 
»” 
4 e (X — X)* — B(Y — V1)? — oZ4** 


where TJ is a term which arises from the limits of integration ac- 


T (3) 


cording to Liebnitz’s rule. 
Now make use of the irrotationality condition that 0w/OX = 
du /d0Z and evaluate as Z —~ 0. Then Eq. (3) becomes 


lim a%w/U) | ~ "aG,(M, Vi) dX, a; 
oe eT an (f) 
Z—-0 ox 4x , (X — Xi) — 8X — ¥,)3}34 


where r now becomes the region on the wing subtended by the 
fore Mach cone of the point Y, Y, 0. 
T—-OasZ —0. 
The boundary condition is 
lim a _" . 
Z 0 (U/Y) = alk, Y) (5) 


Differentiating Eq. (5) with respect to X, 
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lim , : s . 
P er" [O(w/U)/OX] = da/dX (6) 


The order of differentiating and taking the limit has been inter 
changed, which is permissible provided w 
Applying Eq. (6) to Eq. (4), there is obtained, finally, 


“Uh iw 


has a derivative. 


Oa 
ox 


AC,(X1, ¥idX dy 
— X,)? — Bx Y — Y,)?|*” 
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This equation is the basic one rather than Eq. (5) of the paper, 
One can now proceed with the derivation of Eq. (10) by modifying 
each of Eqs. (6) through (10) as follows: 

lim oO : 
2-002 J_« 
each of the double integrals, and substitute 0a/Ox for a. 


Drop the operation dx, which appears in front of 


The 
derivations of Eqs. (34) and (43) should be modified in a similag 
mariner 





Solution of the One-Dimensional Flow Equations of a Viscous, Heat- 
Conducting, Compressible Gas 


(Continued from page 684) 
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Mixing of Any Number of Streams in a Duct of Constant Cross-Sectional 
Area 


(Continued from page 698) 


It is interesting to note that the same final Mach 
Number is obtained with any combination of initial 


area ratios and total pressures that satisfy the above 
specified condition in regions (1) and (2). 
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